knitr::opts_chunk$set(echo = TRUE)
library(data.table)
library(imager)
library(visDec)
library(ggplot2)
library(rpart)
library(rattle)
library(ROCR)
library(caret)
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.

load("results/deBiltResults2015.RData")
ggplot(imageSummary, aes(x = log(MOR), y = meanEdge)) + geom_point() 
imClear <-  subim(load.image(filenames[75]), y > 16)
#```r
#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
#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.

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)
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).

source("scripts/ReadMeteoData.R")
windDeBilt <- ReadWindData("inst/extdata/Sensor/DeBiltWind1-1-2015-31-08-2016.csv")
#windTwente <- ReadWindData()
setkey(windDeBilt, dateTime)
windData <- windDeBilt
humidityDeBilt <- ReadHumidityData("inst/extdata/Sensor/DeBilt_RelativeHumidity.csv")
setkey(humidityDeBilt, dateTime)
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)
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)
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:

fogTree <- rpart(foggy ~ meanEdge + changePoint + meanBrightness + smoothness +
                   meanHue + meanSaturation,
                 trainDeBilt,
                 control = rpart.control(cp = 0.015))
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.

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)
}
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.

fogTreeMeteo <- rpart(foggy ~ meanEdge + changePoint + meanBrightness +
                       windSpeed + relHumidity, 
                     trainDeBilt , control = rpart.control(cp = 0.015))
fancyRpartPlot(fogTreeMeteo, sub="")
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.

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.

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:

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)")
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:

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.

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

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).

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)
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)
source("scripts/ReadMeteoData.R")
windTwente <- ReadWindData("inst/extdata/Sensor/TwenteWind1-1-2015-31-08-2016.csv")
setkey(windTwente, dateTime)
humidityTwente <- ReadHumidityData("inst/extdata/Sensor/TwenteHumidity1-1-2015-31-08-2016.csv")
setkey(humidityTwente, dateTime)
tempDewPointTwente <- ReadTempDewPointData("inst/extdata/Sensor/TwenteTemp_DewPoint1-1-2015-31-08-2016.csv")
setkey(tempDewPointTwente, dateTime)
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)
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:

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

The accurancy of the classification is as follows for the :

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")
mlm <- lm(log(MOR) ~ meanEdge + changePoint + meanBrightness, trainTwente) 
round(mlm$coefficients, 3)
predVals <- predict(mlm, trainTwente)
trainTwente[, predMOR := predVals]
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.

#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.

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.

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.

#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.

#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.

fogTreeMeteoTwente <- rpart(foggy ~ meanEdge + changePoint + meanBrightness +
                       windSpeed + relHumidity + airTemperature + dewPoint, 
                     trainTwente , control = rpart.control(cp = 0.015))
fancyRpartPlot(fogTreeMeteoTwente, sub="")
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

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.

inspectLowHumPics<-trainTwente[meanEdge<0.008352802 & relHumidity<45.5 & meanEdge>=0.004425891 & airTemperature>=1.85 & relHumidity <99.5]

Faulty MOR detection on 23-08-2015

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.

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.

fancyRpartPlot(fogTreeRainTwente, sub="")
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.

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.

fancyRpartPlot(fogTreeBinaryRainTwente, sub="")
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.

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.

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) ]
fogTreeNoRain <- rpart(foggy ~ meanEdge + changePoint + meanBrightness +
                                 relHumidity + windSpeed + airTemperature + dewPoint + precipitationIntElec + 
                                 precipitationIntPWS + precipitationDurationElec + precipitationDurationPWS, 
                                 trainTwenteNoRain , control = rpart.control(cp = 0.015))
fancyRpartPlot(fogTreeNoRain, sub="")

The accurancy of the classification is as follows for the :

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

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.

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

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



MartinRoth/fogDec documentation built on May 7, 2019, 3:38 p.m.