ApproachDescription.Rmd

---
title: "Description of Approach"
author: "Martin Roth and Andrea Pagani"
date: "September 29, 2016"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, Libraries, include=FALSE}
library(data.table)
library(imager)
library(visDec)
library(ggplot2)
library(rpart)
library(rattle)
library(ROCR)
library(caret)
```

```{r, ExamplePictures, include=FALSE ,cache=TRUE}
path <- system.file("extdata/Meetterrein", package="visDec")
filenames <- list.files(path,
                        pattern=glob2rx("Meetterrein_201510*.jpg"),
                        full.names=TRUE)
```

## Global feature approach

One problem with the presented landmark discrimination and contrast reduction
approaches is, that they require specific targets to be selected. An approach
to overcome this by focusing on image features is presented in the
following. First, we present the features we looked at so far and then we
present two different ways to use these features.

### Mean edges
Edge detection was used already in the landmark  discrimination approach. 
Instead of focusing on specific targets, we calculate now the mean number of 
edges in a given picture. Under daylight conditions this number can be viewed as
a relative indication of the fogginess. 

```{r, meanEdgeRegression, cache=TRUE, echo=FALSE, warning=FALSE, fig.height=4, fig.width=4}
load("results/deBiltResults2015.RData")
ggplot(imageSummary, aes(x = log(MOR), y = meanEdge)) + geom_point() 
```

```{r, Edges, cache=TRUE, include=FALSE, fig.height=4, fig.width=8.4}
imClear <-  subim(load.image(filenames[75]), y > 16)
#```{r, EdgesPlotClear, echo=FALSE, fig.cap=paste("Edges in a clear situation"), fig.height=4, fig.width=8.4}
#im <- subim(load.image(filenames[75]), y > 16)
#old_par <- par(mfrow=c(1,2))
#plot(im)
#DetectEdges(im) %>% plot
#par(old_par)
#```
#In a clear situation we observe `r round(DetectMeanEdges(im)*100,2)`% edges.
#
#```{r, EdgesPlotFoggy, echo=FALSE, fig.cap=paste("Edges in a foggy situation"), fig.height=4, fig.width=8.4}
#im <- subim(load.image(filenames[49]), y > 16)
#old_par <- par(mfrow=c(1,2))
#plot(im)
#DetectEdges(im) %>% plot
#par(old_par)
#```
imFoggy <- subim(load.image(filenames[49]), y > 16)
old_par <- par(mfrow=c(1,2))
plot(imClear)
plot(imFoggy)
par(old_par)
```

In the clear situation we have a fraction of 
`r round(DetectMeanEdges(imClear)*100,2)`% edges and in in the foggy situation
we observe `r round(DetectMeanEdges(imFoggy)*100,2)`% edges.

Note that in under night conditions fog actually results in relatively more
edges, because of the scattering of light sources by the water particles.  

### Transmission estimate using the Dark Channel prior

He et al. (2008) present an approach for haze removal from (single) images based
on the so-called dark channel prior. The dark channel prior is based on the
observation, "that most patches (small subimages) in a haze-free outdoor image
contain some pixels which have very low intensities in at least one color
channel". The idea was to follow this approach but instead of focusing on the
haze removal, we focus on the haze instead. Therefore, we estimate the
transmission in the local patches with the technique developed by 
He et al. (2008). From the transmission we estimate horizontal averages to 
infer the transition between sky and non-sky regions. In a clear situation this 
transition is relatively sharp. However, in a foggy situation the transition 
is less clear.

```{r, TransmissionPlotClear, cache = TRUE, echo=FALSE, fig.cap="Transmission in a clear situation", fig.height=3, fig.width=8.4, dependson='ExamplePictures'}
im <- subim(load.image(filenames[75]), y > 16)

darkChannel       <- GetDarkChannel(im)
atmosphere        <- GetAtmosphere(im, darkChannel)
transmissionClear <- GetTransmissionEstimate(im, atmosphere)

old_par <- par(mfrow=c(1,3))
plot(im)
as.cimg(transmissionClear) %>% plot
plot(GetHorizAvgTrans(im), xlab="", ylab="")
par(old_par)
```

```{r, TransmissionPlotFoggy, cache=TRUE, echo=FALSE, fig.cap="Transmission in a foggy situation", fig.height=3, fig.width=8.4, dependson='ExamplePictures'}
im <- subim(load.image(filenames[49]), y > 16)

darkChannel       <- GetDarkChannel(im)
atmosphere        <- GetAtmosphere(im, darkChannel)
transmissionFoggy <- GetTransmissionEstimate(im, atmosphere)

old_par <- par(mfrow=c(1,3))
plot(im)
as.cimg(transmissionFoggy) %>% plot()
plot(GetHorizAvgTrans(im), xlab="", ylab="")
par(old_par)
```

The idea is now to use the horizontal averages and compute the changepoint or a
measure of the variation of the curve and use these as features of the image.

It is clear that different features could be used as well. We looked for 
instance also at different color spaces and in the following we use also the
mean brightness of the image as a feature. 

### Classification
Using these (and possibly other) features we want to
develop a classification (and finally a clustering) scheme, to predict whether a
image is foggy or not. As training set we have data for the measurement site at
De Bilt, from June 1 2015 until December 31 2015 (16001 images). 
The training set are the images for the same fixed camera for the period January
1 till June 30 2016 (13883 images).

```{r, WindData, cache=TRUE}
source("scripts/ReadMeteoData.R")
windDeBilt <- ReadWindData("inst/extdata/Sensor/DeBiltWind1-1-2015-31-08-2016.csv")
#windTwente <- ReadWindData()
setkey(windDeBilt, dateTime)
windData <- windDeBilt
```

```{r, HumidityData, cache=TRUE}
humidityDeBilt <- ReadHumidityData("inst/extdata/Sensor/DeBilt_RelativeHumidity.csv")
setkey(humidityDeBilt, dateTime)
```


```{r, dataDeBilt, include=FALSE, cache=TRUE, dependson=c('WindData','HumidityData')}
dataDeBilt <- readRDS("results/ResultsDeBilt2015-2016_3hSun.rds")
offsetBeforeSunrise <- 0 #time in minutes
offsetAfterSunset <- 0
#trainDeBilt <- imageSummary[dateTime > sunriseDateTime - offsetBeforeSunrise * 60 & dateTime < sunsetDateTime + offsetAfterSunset * 60, ]
dataDeBilt <- dataDeBilt[dateTime > sunriseDateTime & dateTime < sunsetDateTime, ]
dataDeBilt[, foggy := MOR < 250]
dataDeBilt[, foggyFactor := as.factor(ifelse(MOR < 250, 'denseFog', 'notspecified'))]
setkey(dataDeBilt, dateTime)
dataDeBilt <- dataDeBilt[windDeBilt, nomatch = 0]
dataDeBilt <- dataDeBilt[humidityDeBilt, nomatch = 0]
dataDeBilt <- na.omit(dataDeBilt)

trainDeBilt <- dataDeBilt[year==2015, .(dateTime, meanEdge, changePoint, smoothness, meanHue,
                meanSaturation, meanBrightness, MOR, foggy, foggyFactor, windSpeed, relHumidity) ]

testDeBilt <- dataDeBilt[year==2016, .(dateTime, meanEdge, changePoint, smoothness, meanHue,
                meanSaturation, meanBrightness, MOR, foggy, foggyFactor, windSpeed, relHumidity) ]

# features <- c("meanEdge", "changePoint", "smoothness", "meanHue", "meanSaturation", "meanBrightness", "windSpeed")
# 
# 
# StandardizeFeatures <- function(dataSet, features) {
#   dataSet[, (features) := lapply(.SD, function(x) {(x - mean(x))/sd(x) }), .SDcols = features]
#   return(dataSet)
# }
# 
# trainDeBilt <- StandardizeFeatures(trainDeBilt, features)


```

```{r, CaretClassifyDeBilt, cache = TRUE, dependson='dataDeBilt'}
set.seed(123)

#inTrain = createDataPartition(y = dataDeBilt$foggy, p = 0.5, list = FALSE)
#str(inTrain)
#
#training <- dataDeBilt[inTrain, ]
#testing  <- dataDeBilt[-inTrain,]

rfModel <- train(foggyFactor ~ meanEdge + changePoint + smoothness + meanBrightness +
                   meanHue + meanSaturation,
                 data = dataDeBilt,
                 metric = "ROC",
                 method = 'ranger',
                 trControl = trainControl(method = "cv", number = 5, verboseIter = TRUE,
                                          classProbs = TRUE))
plot(rfModel)


```

```{r, CaretConfusionDeBilt, cache = TRUE, dependson='CaretClassifyDeBilt'}

rfClassesTest  <- predict(rfModel, newdata = testDeBilt)
rfClassesTrain <- predict(rfModel, newdata = trainDeBilt)

confusionMatrix(data = rfClassesTrain, trainDeBilt$foggyFactor)
confusionMatrix(data = rfClassesTest, testDeBilt$foggyFactor)

```

Based on the training set we obtain the following preliminary decision tree:
```{r, Classification, cache = TRUE, include=FALSE, dependson='dataDeBilt'}
fogTree <- rpart(foggy ~ meanEdge + changePoint + meanBrightness + smoothness +
                   meanHue + meanSaturation,
                 trainDeBilt,
                 control = rpart.control(cp = 0.015))
```

```{r, ClassificationTree, echo=FALSE, cache=TRUE, dependson='Classification'}
fancyRpartPlot(fogTree, sub="")
```

The tree has to be read as follows for an image with `meanEdge = 0.0039` and 
`changePoint = 295`. From the top node the lower right node,
with number 3, is reached because `meanEdge >= 0.0041` is false. 
In the next step the lower left node, with number 6, is reached 
because `changePoint < 296` is true. 
This node number 6, tells us that there were 24 instances in the training
set and from these 0 were foggy, i.e. dense fog with `MOR < 250`.
Therefore, the chance that an image in this class is foggy is set to zero.



```{r, echo=F}
evaluatePrediction <- function(t){
cleanedTest <- t[complete.cases(t)] #there are sitations where the MOR is not verified and is given an NA for the measure, thus remove, actually we should remove from the train set and test set before giving to the classifier

predBin <- ifelse(cleanedTest$pred > 0.01, 1, 0) #round(cleanedTest$pred,0)

labels <- cleanedTest$foggy
tp <- sum(predBin == 1 & labels == T)
tp
tn <- sum(predBin == 0 & labels == F)
tn
falsePos <- sum(predBin == 1 & labels ==F)
falsePos
falseNeg <- sum(predBin == 0 & labels ==T)
falseNeg
precision <- tp / (tp + falsePos)
recall <- tp / (tp + falseNeg)
print("precision:")
print(precision)
print("recall:")
print(recall)

print("F1 score:")
print(2*precision*recall/(precision+recall))

# predComp<-prediction(cleanedTest$pred, labels)
# PR.perf <- performance(predComp, "acc")
# plot(PR.perf)
}

```

```{r, EvaluatePrediction, cache=TRUE, dependson=c('dataDeBilt', 'Classification')}
hindVals <- predict(fogTree, trainDeBilt, method = "class")
trainDeBilt[, pred := hindVals]
predVals <- predict(fogTree, testDeBilt, method="class")
testDeBilt[, pred := predVals]

confusionMatrix(ifelse(trainDeBilt$pred > 0.3, TRUE, FALSE), trainDeBilt$foggy, positive = "TRUE")
confusionMatrix(ifelse(testDeBilt$pred > 0.3, TRUE, FALSE), testDeBilt$foggy, positive = "TRUE")

```




Other meteorological parameters are added to the set of features analyzed for the machine learning problem.
```{r, MeteoFeatureAddition, cache=TRUE, dependson='dataDeBilt'}

fogTreeMeteo <- rpart(foggy ~ meanEdge + changePoint + meanBrightness +
                       windSpeed + relHumidity, 
                     trainDeBilt , control = rpart.control(cp = 0.015))
```

```{r, ClassTrWithWind, echo=FALSE, cache=TRUE, dependson='MeteoFeatureAddition'}
fancyRpartPlot(fogTreeMeteo, sub="")
```


```{r, TestWithWind, cache=TRUE, dependson=c('dataDeBilt', 'MeteoFeatureAddition')}
hindVals <- predict(fogTreeMeteo, trainDeBilt, method = "class")
trainDeBilt[, pred := hindVals]
predVals <- predict(fogTreeMeteo, testDeBilt, method="class")
testDeBilt[, pred := predVals]

confusionMatrix(ifelse(trainDeBilt$pred > 0.3, TRUE, FALSE), trainDeBilt$foggy, positive = "TRUE")
confusionMatrix(ifelse(testDeBilt$pred > 0.3, TRUE, FALSE), testDeBilt$foggy, positive = "TRUE")

```

In 2015 we have `r trainDeBilt[foggy == TRUE, .N]` images
in the test set, that are labeled as foggy, i.e. a MOR below 250 meters.

We use the obtained tree to predict the probability of fog from the features of
the images in the test set. In 2016 we have `r testDeBilt[foggy == TRUE, .N]` images
in the test set, that are labeled as foggy, i.e. a MOR below 250 meters. The
table below shows the days, where we either had a high probability of fog
(`> 50%`) and a relatively high MOR (`> 250 meters`) or days with a low 
probability of fog (`< 50%`) and a relatively low MOR (`< 250 meters`).
In total we obtain five days for which this condition is true. 

<!--```{r, ClassificationErrors, echo=FALSE}
knitr::kable(cleanedTest[(pred > 0.5 & MOR > 250) | (pred < 0.5 & MOR < 250),.(dateTime, MOR, meanEdge, changePoint, foggy, pred = round(pred, 2))])
-->

<!--
January 1 2016 is a foggy day. The MOR and the classification mostly give the 
same results. However, at 9:30 and 13:10 the change point detection based on the
transmission estimate leads to wrong a conclusion. 
-->
On January 6 it is already
dark (which should have been filtered out prior to the analysis (time zone
issues)). Moreover, the MOR gives values around 750 meters so the roughly 57% 
probability of fog do not seem too bad. 
On January 23 one is not able to see the radar tower, which is certainly closer
than the reported 1.5 km from the MOR, so here the camera based approach gives
the better result. 

The situation in Twente is different for two reasons, the wide angle of the
camera makes the horizontal averaging of the transmission rate less appropriate.
Moreover, even on a clear day there are only a few edges in the image (which are
mostly very close to the camera, i.e. from the equipment of the automatic 
weather station). Nevertheless, it was quite simple to detect failures of the 
visibility sensor using the two described features, such as during the afternoon
of August 23 2015, where the sensor consistently gave MOR < 250, although the 
image is very clear.

### Regression
When we are interested not only in the distinction between foggy and not foggy,
we can try to build a regression model. In a first exploratory approach we model 
`log(MOR)` as a linear function of `meanEdge + changePoint + meanBrightness`.
```{r, MLM, echo=FALSE, cache=TRUE, dependson='dataDeBilt'}
mlm <- lm(log(MOR) ~ meanEdge + changePoint + meanBrightness, trainDeBilt) 
round(mlm$coefficients, 3)
predVals <- predict(mlm, trainDeBilt)
trainDeBilt[, predMOR := predVals]
```

We can plot the modeled response versus the observations, as seen below:

```{r, UncleanPlot, warning=FALSE, echo=FALSE, fig.height=4, fig.width=4, cache=TRUE, dependson='MLM'}
ggplot(trainDeBilt, aes(x = log(MOR), y = predMOR)) + geom_point() +
  geom_vline(xintercept = c(log(250), log(1000), log(3000), log(5000)), lty = 3) + 
  ylab("modelled log(MOR)")
```

```{r, MLMWind, echo=FALSE, cache=TRUE, dependson='dataDeBilt'}
mlm <- lm(log(MOR) ~ meanEdge + changePoint + meanBrightness + windSpeed, trainDeBilt) 
round(mlm$coefficients, 3)
predVals <- predict(mlm, trainDeBilt)
trainDeBilt[, predMOR := predVals]
```

We can plot the modeled response versus the observations, as seen below:

```{r, UncleanPlotWind, warning=FALSE, echo=FALSE, fig.height=4, fig.width=4, cache=TRUE, dependson='MLMWind'}
ggplot(trainDeBilt, aes(x = log(MOR), y = predMOR)) + geom_point() +
  geom_vline(xintercept = c(log(250), log(1000), log(3000), log(5000)), lty = 3) + 
  ylab("modelled log(MOR)")
```

There are some quite apparent features in this plot. For instance the points 
with the very low MOR and a prediction of around 10 are mostly due to an 
failure of the transmissometer. The values with the largest prediction on the
other hand correspond to different sceneries in the picture (we believed the
camera was fixed) or situations where it is already dark. The point on the lower
right corner with a prediction below 7 corresponds again to a different scenery.

```{r, CleanData, echo=FALSE, cache = TRUE, dependson='MLM'}
trainDeBilt[, id := 1:.N]
#
# ggplot(train, aes(x = log(MOR), y = predMOR)) + geom_point()
#
# train[predMOR > 12]
# June 22 2015 11:10 different scenery
# August 4 2015 11:00 different scenery
# September 27 19:10 / 19:20 already dark 
# October 20 18:00 already dark
#
trainDeBilt <- trainDeBilt[predMOR < 12]
#
# ggplot(train, aes(x = log(MOR), y = predMOR)) + geom_point()
#
# train[(log(MOR) < 6 & predMOR > 9)]
# November 1 2015 clear view (afer foogy period / might be time issue)
# November 2 2015 clear view (no fog)
#
trainDeBilt <- trainDeBilt[!(log(MOR) < 6 & predMOR > 9)]
#
# ggplot(train, aes(x = log(MOR), y = predMOR)) + geom_point()
#
# train[(log(MOR) > 8 & predMOR < 9)]
# June 19 2015 9:20 different scenery
# otherwise already dark
#
trainDeBilt <- trainDeBilt[!(log(MOR) > 8 & predMOR < 9)]
#
# ggplot(train, aes(x = log(MOR), y = predMOR)) + geom_point()
#
# train[(log(MOR) < 7 & predMOR > 9.5)]
# June 5 2015 too late
# September 23 too late
# November 2 clear view
# November 26 clear view
#
trainDeBilt <- trainDeBilt[!(log(MOR) < 7 & predMOR > 9.5)]
#
# ggplot(train, aes(x = log(MOR), y = predMOR)) + geom_point()
#
# train[(log(MOR) < 8 & predMOR > 10.2)]
# June too late
# July too late
# November 1 until 16 clear view
# November 1 17:10 already dark
# November 2 29:10 clear view
# November 26 8:40 clear view
#
trainDeBilt <- trainDeBilt[!(log(MOR) < 8 & predMOR > 10.2)]
#
# ggplot(train, aes(x = log(MOR), y = predMOR)) + geom_point()
#
# train[(log(MOR) < 10 & predMOR > 11)]
# June 22 already too late
# October 30 Car in the scenery
#
trainDeBilt <- trainDeBilt[!(log(MOR) < 10 & predMOR > 11)]
#
#
#
```

#### Outlook

From the previous section, it is clear that the underlying pictures have to be
reviewed carefully. In the end we want to be able to allow for different views
of the camera, but for this exploratory analysis we have to exclude these data
from the analysis.

Regarding cameras at sites, where no MOR sensor is available regression seems a
big challenge. But the idea of clustering the data, maybe in a semi supervised
way (as there are usually only a few foggy situations), and reporting if the 
probability of fog exceeds a certain level seems relatively close.

#### Outlier detection
```{r Outliers, dependson='dataDeBilt', eval=FALSE}
dataDeBilt <- dataDeBilt[changePoint < 400]
dataDeBilt <- dataDeBilt[changePoint > 167]
```

### Twente test field

The images from the test field in Twente present some additional challenges 
when to be elaborated.
First of all the images are taken from a "Fish-eye" camera (A fisheye lens is 
an ultra wide-angle lens that produces strong visual distortion intended to 
create a wide panoramic or hemispherical image. Fisheye lenses achieve extremely
wide angles of view by forgoing producing images with straight lines of 
perspective (rectilinear images), opting instead for a special mapping 
(for example: equisolid angle), which gives images a characteristic convex 
non-rectilinear appearance. Source: wikipedia).

After a preliminary analysis we noted that the metrics characterizing the image
(i.e., mean edges, changepoint) were also influenced by the wide angle 
objective. We have then proceeded to rectify the images accordingly. 
The rectification has been done by using an ad-hoc script leveraging on the 
internal capabilities of GNU Image Manipulation Program (GIMP).

```{r, ImagesExTw, cache=TRUE, echo=FALSE}
exampleImg <- c("inst/extdata/EHTW_201508011440.jpg","inst/extdata/EHTW_201508011440_Rectified.jpg")
example_par <- par(mfrow=c(1,2))
img <- lapply(exampleImg, load.image)
plot(img[[1]], main = "Original Fish-eye image")
plot(img[[2]], main = "Rectified image")
par(example_par)
```


```{r, dataTw, cache=TRUE}
dataTwente <- readRDS("results/ResultsTwente2015-2016_3hSun.rds")
offsetBeforeSunrise <- 0 #time in minutes
offsetAfterSunset <- 0
#trainDeBilt <- imageSummary[dateTime > sunriseDateTime - offsetBeforeSunrise * 60 & dateTime < sunsetDateTime + offsetAfterSunset * 60, ]
dataTwente <- dataTwente[dateTime > sunriseDateTime & dateTime < sunsetDateTime, ]
dataTwente[, foggy := MOR < 250]
setkey(dataTwente, dateTime)
```

```{r, WindTw, cache=TRUE}
source("scripts/ReadMeteoData.R")
windTwente <- ReadWindData("inst/extdata/Sensor/TwenteWind1-1-2015-31-08-2016.csv")
setkey(windTwente, dateTime)
```

```{r, HumidityTw, cache=TRUE}
humidityTwente <- ReadHumidityData("inst/extdata/Sensor/TwenteHumidity1-1-2015-31-08-2016.csv")
setkey(humidityTwente, dateTime)
```

```{r, TempTw, cache=TRUE}
tempDewPointTwente <- ReadTempDewPointData("inst/extdata/Sensor/TwenteTemp_DewPoint1-1-2015-31-08-2016.csv")
setkey(tempDewPointTwente, dateTime)
```

```{r, PrecipTw, cache=TRUE}
precipitationTwente <- ReadPrecipitationData("inst/extdata/Sensor/TwenteRainAll-1-2015-31-08-2016.csv")
precipitationTwente[, rain := FALSE]
rainDuration <- 400
precipitationTwente[ precipitationDurationPWS > rainDuration | precipitationDurationElec > rainDuration, rain:= TRUE]
setkey(precipitationTwente, dateTime)
```



```{r, trainTestTw, cache=TRUE, dependson=c('TempTw', 'HumidityTw', 'WindTw', 'dataTw', 'PrecipTw')}
dataTwente <- SynchronizeSensorReadingsNoMORPicture(windTwente, dataTwente)
dataTwente <- SynchronizeSensorReadingsNoMORPicture(humidityTwente, dataTwente)
dataTwente <- SynchronizeSensorReadingsNoMORPicture(tempDewPointTwente, dataTwente)
dataTwente <- SynchronizeSensorReadingsNoMORPicture(precipitationTwente, dataTwente)


dataTwente <- na.omit(dataTwente)

trainTwente <- dataTwente[year==2015, .(dateTime, meanEdge, changePoint, smoothness, meanHue,
                meanSaturation, meanBrightness, MOR, foggy, day, month, hour, windSpeed, relHumidity, airTemperature, dewPoint,
                precipitationIntElec, precipitationIntPWS, precipitationDurationElec, precipitationDurationPWS, rain) ]

testTwente <- dataTwente[year==2016, .(dateTime, meanEdge, changePoint, smoothness, meanHue,
                meanSaturation, meanBrightness, MOR, foggy, day, month, hour, windSpeed, relHumidity, airTemperature, dewPoint,
                precipitationIntElec, precipitationIntPWS, precipitationDurationElec, precipitationDurationPWS, rain) ]

```

In 2015 we have `r trainTwente[foggy == TRUE, .N]` images
in the test set, that are labeled as foggy, i.e. a MOR below 250 meters.

Based on the training set we obtain the following preliminary decision tree 
using just the features of the images:
```{r, ClassTw, include=FALSE, cache=TRUE, dependson='trainTestTw'}


fogTree <- rpart(foggy ~ meanEdge + changePoint + meanBrightness, trainTwente , control = rpart.control(cp = 0.019))
```

```{r, ClassTreeTw, echo=FALSE, cache=TRUE, dependson='ClassTw'}
fancyRpartPlot(fogTree, sub="")
```

The accurancy of the classification is as follows for the :
```{r, EvalPredTw, cache=TRUE, dependson='ClassTw'}
hindVals <- predict(fogTree, trainTwente, method = "class")
trainTwente[, pred := hindVals]

predVals <- predict(fogTree, testTwente, method="class")
testTwente[, pred := predVals]

#usually it is not a good practice to show the results of the classification on the training set itself.
confusionMatrix(ifelse(trainTwente$pred > 0.3, TRUE, FALSE), trainTwente$foggy, positive = "TRUE")
confusionMatrix(ifelse(testTwente$pred  > 0.3, TRUE, FALSE), testTwente$foggy, positive = "TRUE")

```


```{r, MLRegTw, echo=FALSE, cache=TRUE, dependson='EvalPredTw'}
mlm <- lm(log(MOR) ~ meanEdge + changePoint + meanBrightness, trainTwente) 
round(mlm$coefficients, 3)
predVals <- predict(mlm, trainTwente)
trainTwente[, predMOR := predVals]
```

```{r, UncleanPlotTwente, warning=FALSE, echo=FALSE, fig.height=4, fig.width=4, cache=TRUE, depenson='MLRegTw'}
ggplot(trainTwente, aes(x = log(MOR), y = predMOR)) + geom_point() +
  geom_vline(xintercept = c(log(250), log(1000), log(3000), log(5000)), lty = 3) + 
  ylab("modelled log(MOR)")
```



###Corner cases situations

There are several situations that pose challenges in the correct classification
of fog conditions and the estimation of the visibility range. In Twente, the
situation is even more complex and more difficult for the classifier given
the environmental conditions.  

```{r, LookIntoClusters, cache=TRUE, echo=FALSE}
#predwrong <- trainTwente[predMOR >9 & log(MOR)<6]
# May 24 2015 04:20 surface fog in an previous time, picture looks clear
# April 6 2015 03:50 slightly surface fog
# August 23 14:00 / 18:40 sensor failure, clear day 
# September 9 05:10 attenuation of fog (horizon visible) in a denser foggy 
#period before and after
#September 26 05:30 direct ray of light in the objective for a few shots 
#might modify image properties
#November 02 10:50 fog starts to dissolve (delay  in taking picture and read 
#of MOR in KIS?)

#predwrong2 <- trainTwente[predMOR >11 & log(MOR)<9.5]
#March 31 2015 16:35-45 dark light and colors (jammed colors?) in the image 
#very different from other images previous next
#April 28 2015 05:15-35 overexposed images light into objective
```
  
First, the objective is unprotected against direct light which has important 
consequences for distinction of objects in the background, the appearance of 
fake objects due to light spots.  

```{r, LookIntoClusters1, cache=TRUE, echo=FALSE}
exampleImg <- c("inst/extdata/cornerCasesTwente/EHTW_201503240815.jpg")
example_par <- par(mfrow=c(1,2))
image <- load.image(exampleImg)
plot(image, main = "Light in the objective")
#plot(0,0)
par(example_par)
```
  
Second, it seems that some humidity is sometimes trapped into the objective 
which gives rise to lines of different colors in the image through the whole 
day until evaporation.  


```{r, LookIntoClusters2, cache=TRUE, echo=FALSE}
exampleImg2 <- c("inst/extdata/cornerCasesTwente/EHTW_201510271350.jpg")
example_par2 <- par(mfrow=c(1,2))
img3 <- load.image(exampleImg2)
plot(img3, main = "Umidity lines on the objective", cex.main = 0.7)
#plot(0,0)
par(example_par2)
```

  
Third, there are objects that appear in the scene and that might cause problems
in the recognition of edges. A dark brown land line appears at some time maybe 
due to some deforestation or enablings works of the countryside; sometimes 
distant vehicles appear on the horizon; wooden poles are installed at the 
perimeter of the weather station.

  
```{r, LookIntoClusters3, cache=TRUE, echo=FALSE}

#predwrong3 <- trainTwente[predMOR >10 & log(MOR)<8.5 & log(MOR)>8]
#March 24 2015 08:05-09:05 light direct to the objective horizon too bright to be recognized
#March 28 2015 05:25-05:55 some light into the objective (colors not sharp)
#April 07 2015 05:55-06:05 some light into the objective (colors not sharp) as above
#April 09 2015 05:55-06:15 light direct to the objective horizon to bright to be recognized left part of image
#April 19 2015 05:15 OK, colors not so sharp
#April 20 2015 05:25-05:45 light direct to the objective horizon to bright to be recognized left part of image
#April 21 2015 06:25 light direct to the objective horizon to bright to be recognized left part of image
#April 24 2015 04:45-06:05 light direct to the objective horizon to bright to be recognized left part of image
#October 27 2015 10:50-14:30 some water drops paths seems to be present on the camera



exampleImg <- c("inst/extdata/cornerCasesTwente/EHTW_201504120955.jpg","inst/extdata/cornerCasesTwente/EHTW_201508251300.jpg")
example_par <- par(mfrow=c(1,2))
img2 <- load.image(exampleImg[[1]])
img3 <- load.image(exampleImg[[2]])
plot(img2, main = "Change in scene: light brown pattern in the distance", cex.main = 0.7)
plot(img3, main = "Change in scene: full green color in the distance", cex.main = 0.7)
par(example_par)


exampleImg2 <- c("inst/extdata/cornerCasesTwente/EHTW_201509170520.jpg")
example_par <- par(mfrow=c(1,2))
img3 <- load.image(exampleImg2)
plot(img3, main = "Objects appearing in the scene: working vehicles in the distance", cex.main = 0.7)
par(example_par)


exampleImg <- c("inst/extdata/cornerCasesTwente/EHTW_201512211403.jpg")
example_par <- par(mfrow=c(1,2))
img2 <- load.image(exampleImg)
plot(img2, main = "Change in scene: wooden poles", cex.main = 0.7)
#plot(0,0)
par(example_par)
```

  
Other more complex conditions that almost always bring to False Positive 
classification are due to the presence of high humidity in the objective or 
drops of rains on the objective protection lid. The camera has unfortunately 
no protection against the rain.

  
  
```{r, LookIntoClusters4, cache=TRUE, echo=FALSE}

#predwrong4 <- trainTwente[predMOR <10 & log(MOR)>10.75]
#March 21 2015 16:25-17:45 quite dark colors and light
#March 22 2015 05:25-05:55 quite dark colors and light
#March 26 2015 06:25-06:55 Ok, but might be quite dark colors and light
#April 12 2015 17:25-18:15 Ok, scene change compared to normal appearance of a light brown piece of land after the woods (cut of bushes/grass?). In certain light conditions it is more evident
#April 16 2015 16:35-18:25 Ok, scene change compared to normal appearance of a light brown piece of land after the woods (cut of bushes/grass?). In certain light conditions it is more evident
#June 20 2015 12:00-19:50 Ok
#July 21 2015 19:20 long shadows from sun in the back
#August 24 2015 14:10-17:10 Ok
#September 17 2015 05:10-07:20 some objects are appearing in the back in the scenes (vehicles working on the land)
#December 22 2015 13:23-15:03 Ok, wooden poles around the measurement field are present in the scene, but were already present since october

#predwrong5 <- trainTwente[predMOR >10 & log(MOR)>6 & log(MOR)<7]
#October 26 2015 13:40-15:10 some shadows on the objective (due to previous drops?)
#August 23 2015 13:00-13:50 clear day (will follow a sensor failure)


exampleImg2 <- list.files("inst/extdata/cornerCasesTwente/", pattern = "EHTW_201506280400.jpg", full.names = TRUE)
example_par <- par(mfrow=c(1,2))
img3 <- load.image(exampleImg2)
plot(img3, main = "Humidity in the objective", cex.main = 0.7)
par(example_par)



exampleImg2 <- c("inst/extdata/cornerCasesTwente/EHTW_201510060700.jpg")
example_par <- par(mfrow=c(1,2))
img3 <- load.image(exampleImg2)
plot(img3, main = "Water drops on the objective", cex.main = 0.7)
par(example_par)





#predwrong6 <- trainTwente[predMOR >8  & predMOR <9 & log(MOR)<9]
#March 20 2015 05:45-06:55 mist
#June 28 2015 03:20-04:10 FALSE POSITIVE condition due to humidity in the objective image is fully blurred
#July 20 2015 03:40-04:10 FALSE POSITIVE condition due to humidity in the objective image is fully blurred
#August 05 2015 04:00-04:40 FALSE POSITIVE and one detection condition due to humidity in the objective image is fully blurred
#August 29 2015 04:40-05:10 FALSE POSITIVE condition due to humidity in the objective image is fully blurred
#September 09 2015 05:30-06:40 mist
#September 13 2015 15:00 water drops on the objective
#September 16 2015 08:20 FALSE POSITIVE drops of water on the objective image is fully blurred
#September 23 2015 06:20-06:40 mist 1 FALSE POSITIVE, but it is actually misty (260m visib)
#October 4 2015 05:40-06:30 mist conditions 2 FALSE POSITIVE due also to humidity in the camera objective
#October 6 2015 05:50-08:10 water on the objective 3 FALSE POSITIVE
#October 19 2015 08:30-08:50 mist
#October 20 2015 06:10-07:00 mist
#November 02 2015 15:30 mist FALSE POSITIVE it's misty (MOR reports ~300m)
#November 10 2015 09:40 water on the camera
#November 15 2015 13:10 water on the camera
#December 16 2015 water on the camera almost whole day

#predwrong7 <- trainTwente[predMOR <9 & log(MOR)>9]
#September 16 2015 06:30 07:10 07:40 08:30 drops of water on the objective image is fully blurred FALSE POSITIVE (07:40)
#September 22 2015 13:30-13:40 drops of water on the objective image is fully blurred
```

  
Other meteorological parameters are added to the set of features analyzed for 
the machine learning problem for Twente.

  
  
  
```{r, MeteoFeatTw, cache=TRUE, dependson='trainTestTw'}

fogTreeMeteoTwente <- rpart(foggy ~ meanEdge + changePoint + meanBrightness +
                       windSpeed + relHumidity + airTemperature + dewPoint, 
                     trainTwente , control = rpart.control(cp = 0.015))
```

```{r, ClassTrWindTw, echo=FALSE, cache=TRUE, dependson='MeteoFeatTw'}
fancyRpartPlot(fogTreeMeteoTwente, sub="")
```


```{r, TestWithWindTw, cache=TRUE, dependson=c('trainTestTw', 'MeteoFeatTw')}
hindVals <- predict(fogTreeMeteoTwente, trainTwente, method = "class")
trainTwente[, pred := hindVals]
predVals <- predict(fogTreeMeteoTwente, testTwente, method="class")
testTwente[, pred := predVals]

confusionMatrix(ifelse(trainTwente$pred > 0.3, TRUE, FALSE), trainTwente$foggy, positive = "TRUE")
confusionMatrix(ifelse(testTwente$pred > 0.3, TRUE, FALSE), testTwente$foggy, positive = "TRUE")

```

The results of the test set for Twente show poor values of precision and recall
```{r, cache=TRUE, dependson='TestWithWindTw'}
FP<-testTwente[pred>0.3 & foggy==FALSE]

datesFP <- unique(format(FP$dateTime,"%Y-%m-%d"))


FN<-testTwente[pred<0.3 & foggy==TRUE]

datesFN <- unique(format(FN$dateTime,"%Y-%m-%d"))

```
False positive conditions occur in the following dates `r nrow(datesFP)`

2016-01-04: drops of rain are in the objective of the camera, pictures are 
blurred. Affected pictures 28  
2016-01-05: there is humidity inside the objective. Likely from the previous 
day rain. Affected images 28  
2016-01-06: there is mist in the distance on of the picture and some pictures 
have water drops on the objective. Affected images 25  
2016-01-07: there are sort of zones with stripes of different lighter color, 
it seems again humidity in the objective enclosure. In conditions of low light 
is more evident a fog-like prediction. Images affected 4.  
2016-01-15: there are sort of zones with stripes of different lighter color, 
it seems again humidity in the objective enclosure. In conditions of low light 
is more evident a fog-like prediction. Images affected 4.  
2016-01-17: snow on the ground, less contrast/edges. Images affected 2.  
2016-02-14: water and snow on the objective. Images affected 18.  
2016-02-19: there are sort of zones with stripes of different lighter color, 
it seems again humidity in the objective enclosure. In conditions of low light 
is more evident a fog-like prediction. Images affected 4.  
2016-03-04: water and snow on the objective. Images affected 4.  
2016-03-05: misty conditions. MOR is slightly above 250m. Images affected 5.  
2016-03-11: misty conditions. MOR is slightly above 250m. Images affected 1.  
2016-04-09: conditions are fine, no mist. Images affected 15.  
2016-04-10: conditions are fine, no mist. Images affected 5.  
2016-04-12: conditions are fine, no mist. Images affected 6.  
2016-04-14: conditions are fine, no mist. Images affected 1.  

False negative conditions occur in the following dates  `r nrow(datesFN)`
2016-02-26: conditions are foggy and snow/ice on the ground. Images affected 6.  
2016-03-05: mist. Images affected 1.  
2016-03-15: mist. Images affected 5.  
2016-04-13: there is a ray of light coming straight into the objective. Images 
affected 13.  
2016-05-23: fog only in the surface. Images affected 1.  
2016-06-04: fog condition dissolving. Images affected 8. (1 image the horizon 
is well visible)
2016-06-28: camera is not fully sharp, but fog seems to have dissolved. Images 
affected 4.  

  

We investigated with more attention the result of the decision tree were the 
classification of the training set suggests a high chance of fog in the 
following condition:
  
meanEdge<0.008352802 & relHumidity<45.5 & meanEdge>=0.004425891 & 
airTemperature>=1.85 & relHumidity <99.5
  
Actually after a visual inspection of the 9 of the 10 selected that the MOR 
sensor data reported as foggy, we conclude that the MOR sensor was in a faulty
condition since it is a clear afternoon day as shown in the figures below. 

```{r, cache=TRUE, dependson='TestWithWindTw'}
inspectLowHumPics<-trainTwente[meanEdge<0.008352802 & relHumidity<45.5 & meanEdge>=0.004425891 & airTemperature>=1.85 & relHumidity <99.5]
```

Faulty MOR detection on 23-08-2015  
```{r, cache=TRUE, echo=FALSE}
fautlyIm <- c("inst/extdata/MORError/EHTW_201508231700.jpg")
parFaulty <- par(mfrow=c(1,2))
img <- load.image(fautlyIm)
plot(img)
par(parFaulty)
```



Rain on the camera objective is an issue that creates false positive situations
in Twente. The idea is then to try to add precipitation variables to the 
features to train the decision tree in order to identify and classify 
these situations automatically.
```{r, PrecTrainTrTw, cache=TRUE, echo=FALSE, dependson='trainTestTw'}


fogTreeRainTwente <- rpart(foggy ~ meanEdge + changePoint + meanBrightness +
                                 relHumidity + windSpeed + airTemperature + dewPoint + precipitationIntElec + 
                                 precipitationIntPWS + precipitationDurationElec + precipitationDurationPWS, 
                                 trainTwente , control = rpart.control(cp = 0.015))

```


No feature concerning rain is actually selected in the decision tree.
```{r, ClassTrWithPrecTw, echo=FALSE, cache=TRUE, dependson='PrecTrainTrTw'}
fancyRpartPlot(fogTreeRainTwente, sub="")
```


```{r, ResPrec, dependson=c('PrecTrainTrTw', 'trainTestTw')}
hindVals <- predict(fogTreeRainTwente, trainTwente, method = "class")
  trainTwente[, pred := hindVals]
  predVals <- predict(fogTreeRainTwente, testTwente, method="class")
  testTwente[, pred := predVals]
confusionMatrix(ifelse(trainTwente$pred > 0.3, TRUE, FALSE), trainTwente$foggy, positive = "TRUE")
confusionMatrix(ifelse(testTwente$pred > 0.3, TRUE, FALSE), testTwente$foggy, positive = "TRUE")
```

Therefore a possibility is to create a binary variable that assesses the 
presence of rain instead of focusing  on precise measure of rain amount an 
duration. The binary rain variable is created by looking at the rain duration. 
If one of the sensors measuring the rain duration measures a rain duration of 
more than `r rainDuration` seconds then the binary rain variable is set to TRUE.



```{r, BinPrecTrainTr, cache=TRUE, echo=FALSE, dependson='trainTestTw'}


fogTreeBinaryRainTwente <- rpart(foggy ~ meanEdge + changePoint + meanBrightness +
                                 relHumidity + windSpeed + airTemperature + dewPoint + rain, 
                                 trainTwente , control = rpart.control(cp = 0.015))

```

Also in this case the binary variable is not selected in the decision tree.
```{r, ClassTrWithBinPrecTw, echo=FALSE, cache=TRUE, dependson='BinPrecTrainTr'}
fancyRpartPlot(fogTreeBinaryRainTwente, sub="")
```


```{r, ResBinPrec, echo=FALSE, cache=TRUE, dependson=c('BinPrecTrainTr', 'trainTestTw')}
hindVals <- predict(fogTreeBinaryRainTwente, trainTwente, method = "class")
  trainTwente[, pred := hindVals]
  predVals <- predict(fogTreeBinaryRainTwente, testTwente, method="class")
  testTwente[, pred := predVals]
confusionMatrix(ifelse(trainTwente$pred > 0.3, TRUE, FALSE), trainTwente$foggy, positive = "TRUE")
confusionMatrix(ifelse(testTwente$pred > 0.3, TRUE, FALSE), testTwente$foggy, positive = "TRUE")
```

Inspecting the false positives caused by water drops in the  objective  we can 
recognize that this situation takes place when there is a sustained period of 
rain (it is reported rainy for a full 10-minute period in a 10 minutes sample) 
and the wind is higher than 3.5 m/s.
We remove then such conditions from the dataset and we also remove the day when
the MOR was not working properly giving visibility below 250m but visual 
inspection showed perfect conditions.
```{r, cleanRainTw, echo=FALSE, cache=TRUE, dependson='trainTestTw'}
noRainDataTwente <- dataTwente[precipitationDurationElec<600 & windSpeed<3.5,]

#noRainDataTwente <- noRainDataTwente[day!=23 & month!=08 & year!=2015,]

removedFoggyRainOrHighWind <- dataTwente[foggy==TRUE, .N] - noRainDataTwente[foggy==TRUE, .N]
datesOfFoggyRainOrHighWind <- dataTwente[dateTime%in%setdiff(dataTwente[foggy==TRUE]$dateTime,noRainDataTwente[foggy==TRUE]$dateTime)]$dateTime
```

The dataset is now composed of `r noRainDataTwente[,.N]` pictures. 
With such a filtering procedure also `r removedFoggyRainOrHighWind` foggy 
pictures have been removed from the dataset. 
The pictures removed that are considered foggy in conditions of rain or high 
wind are from `r datesOfFoggyRainOrHighWind`. 
This timestamp is the same as the faulty MOR readings reported above.


We partition again the the dataset and train a decision tree with image and 
meteorological features.

```{r, trainTestNoRainTw, echo=FALSE, cache=TRUE, dependson='cleanRainTw'}

trainTwenteNoRain <- noRainDataTwente[year==2015, .(dateTime, meanEdge, changePoint, smoothness, meanHue,
                meanSaturation, meanBrightness, MOR, foggy, day, month, hour, windSpeed, relHumidity, airTemperature, dewPoint,
                precipitationIntElec, precipitationIntPWS, precipitationDurationElec, precipitationDurationPWS, rain) ]

testTwenteNoRain <- noRainDataTwente[year==2016, .(dateTime, meanEdge, changePoint, smoothness, meanHue,
                meanSaturation, meanBrightness, MOR, foggy, day, month, hour, windSpeed, relHumidity, airTemperature, dewPoint,
                precipitationIntElec, precipitationIntPWS, precipitationDurationElec, precipitationDurationPWS, rain) ]

```




```{r, ClassTwNoRain, include=FALSE, cache=TRUE, dependson='trainTestNoRainTw'}


fogTreeNoRain <- rpart(foggy ~ meanEdge + changePoint + meanBrightness +
                                 relHumidity + windSpeed + airTemperature + dewPoint + precipitationIntElec + 
                                 precipitationIntPWS + precipitationDurationElec + precipitationDurationPWS, 
                                 trainTwenteNoRain , control = rpart.control(cp = 0.015))
```

```{r, ClassTrTwNoRain, echo=FALSE, cache=TRUE, dependson='ClassTwNoRain'}
fancyRpartPlot(fogTreeNoRain, sub="")
```

The accurancy of the classification is as follows for the :
```{r, EvalPredTwNoRain, cache=TRUE, dependson=c('ClassTwNoRain','trainTestNoRainTw')}
hindVals <- predict(fogTreeNoRain, trainTwenteNoRain, method = "class")
trainTwenteNoRain[, pred := hindVals]

predVals <- predict(fogTreeNoRain, testTwenteNoRain, method="class")
testTwenteNoRain[, pred := predVals]

#usually it is not a good practice to show the results of the classification on the training set itself.
confusionMatrix(ifelse(trainTwenteNoRain$pred > 0.3, TRUE, FALSE), trainTwenteNoRain$foggy, positive = "TRUE")
confusionMatrix(ifelse(testTwenteNoRain$pred  > 0.3, TRUE, FALSE), testTwenteNoRain$foggy, positive = "TRUE")

```

This additional filtering helps partially in the reduction of the FALSE 
POSITIVE cases. With a threshold to discriminate between fog and non-fog 
conditions equal 30%, the classifier achives  on the test set a balanced 
accurancy of 81% compared when the rain situations are filtered out, while 
the balanced accuracy was 62% on the non-filtered dataset. Some FALSE POSITIVE 
cases still remain and they are still in situations where drops of rain or 
humidity are present in the objective.


From a visual inspection of the cases with the highest

```{r, cache=TRUE, dependson='trainTestNoRainTw'}
inspectPicts1<-trainTwenteNoRain[meanEdge>=0.006225944 & airTemperature< 2.45 & meanEdge< 0.007368848]

inspectPicts1$dateTime
```
20-03-2015: dissolvng fog
21-04-2015: is this fog or UFO?? strong light in the objective
02-05-2015: light surface fog
02-10-2015: UFOs are back? strong light in the objective
03-11-2015: dissolving fog.

All these situation there is quite a lot of kind-of reflected light.

```{r, cache=TRUE, dependson='trainTestNoRainTw'}
inspectPicts2<-trainTwenteNoRain[meanEdge<0.006225944 & airTemperature>=6.25 & meanBrightness>=0.0020341]

inspectPicts2$dateTime
```
09-09-2015: fog dissolving
02-11-2015: fog forming
19-10-2015: fog
06-10-2015: drop of water

```{r, cache=TRUE, dependson='trainTestNoRainTw'}
inspectPicts3<-trainTwenteNoRain[meanEdge<0.006225944 & airTemperature>=6.25 & meanBrightness<0.0020341]

inspectPicts3$dateTime
```

23-06-2015: humidity in the camera
03-07-2015: forming fog
20-07-2015: humidity/rain in the objective
23-09-2015: foggy sometimes and humidity in objective
27-10-2015: kind of foggy
05-08-2015: lot of humidity in the objective
31-08-2015: high bright sky in the picture involved