Category Archives: Programming

Programming is my life! As I’m working as a programmer, it more or less is. And I’m also coding in my free-time.
Languages I’m proficient in include R, Java, Python and C++. But I’m also capable in writing Bash or Perl scripts, even though I don’t like the second one that much.

In this category of my blog I will generally present code to you. Be it code in form of the implementation of an algorithm or an application. The two languages Iyou will probably find most of here are R and C++. Also because they are easy to integrate into each other. But this might change over time, as programming languages and their use change.

If you want some algorithm explained, just write me. I will see, what I can do!

My Humble Beginnings As A Programmer


Warning: count(): Parameter must be an array or an object that implements Countable in /homepages/10/d771928198/htdocs/clickandbuilds/RandomThoughtsScienceandProgramming/wp-content/plugins/kblog-metadata/kblog-author.php on line 332

Often when I tell people I haven’t really started programming till my time at the university, they won’t believe me. And while I performed well at programming tasks from the beginning of my life as a Bioinformatics student, this is mostly true. But I’m not here today to talk about that, I’m here to talk about the other half of the truth. The programming I did before I was a student or in other words my humble beginnings as a programmer

The Schooldays

Well. As I’m a little bit older, there hasn’t been much programming during my schooldays. Actually, the only time we programmed sometime was in ninth grade in our math class. We had an old teacher, who was an early computer enthusiast. And so when we were done with the subject matter for the year, we did some programming in Turbo Pascal. If you don’t know, Turbo Pascal is an out of fashion programming language, that as far as I know was developed for teaching purposes.

Although we didn’t do great things in it, we at least learnt some concepts. As far as I remember we just had to implement a program, that printed out all the primes from one to one hundred. Because I was done with it pretty quickly, I brought some programming books from my dad to class. Those books contained code for computer games written in Turbo Pascal. Yep… That was actually a thing in the 80’s to 90’s! But I didn’t get far. Do you know how bothersome it is to copy paste pages of code the old school way?
So yep… That wasn’t really my beginnings as a programmer. I mean, I learnt something, but it didn’t really get me motivated to continue learning on my own. And a little bit of time passed till my next attempt.

Beginnings As A Programmer On My Own

When I was about eighteen, the story continued. I was just about finishing middle school, having my first girlfriend and so on. Besides that I played the MMORPG World of Warcraft. But while I started playing on the official servers, at some point of time I started playing on a private shard called Lordearon. For the Warcraft fans amongst you: NO that was not a writing error at my cost. That was actually the name of the server. Playing there it didn’t take long for me to become a gamemaster (GM). For your information, Lordearon was run as a project of some students. So they didn’t have much time and resources and therefore needed help in managing the server.

During this time I learnt a little bit of PHP, which was used in conjunction with SQL to run the server. And to fix bugs it was useful to able to write small scripts. I didn’t really do this for long, but still I learnt some stuff mostly on my own. And that some computer science students trusted me enough to give me GM and even some admin rights, gave me self confidence.

Self developed picture of me, when I was nineteen.

Just Playing Around

After some time had passed I got the idea, that I had to learn C. Well, it was shortly after I broke up with my then-girlfriend. So I had a lot of newly gained free time and this wasn’t my only project I started during this time. I found a nice online tutorial for learning C back then, which I followed a bit. And then I programmed some easy programs for doing some calculations for the laboratory… Nothing special, just simple arithmetic.

Unfortunately I didn’t follow those beginnings as a programmer furthermore. But even if I did only the basics and no algorithms and stuff, I think I laid down a valuable foundation. Then came my second-chance education, where I did my high school degree. During this time I was used to capacity with school. However I improved my math capabilities a lot along the way, which is important for programming.

Before My Student Life

After I finished my high school degree (Abitur) I had some months of free time until my time at the university began. And luckily I choose to spend some of it on online courses. While some of them were about statistics and maths in general, others were about programming basics. And I noticed that I should’ve started with something like that much earlier, because I liked it. Again, I didn’t learn any hard concepts, but those basics helped me to focus on the hard problems later at the university.

So this is also, what I would advise you, if you’re planning on studying a computer science related program. You don’t need to know the hard concepts, but most lectures won’t give you the time to understand the basics of programming like loops and conditional statements. Additionally, it is also of advantage to know how to read and write data. That’s something I didn’t learn before uni. While those are rather easy concepts, if you compare them to other ones you learn as a computer scientist, if you have never done such a thing, you could struggle a bit at first.

Nowadays there’s really this big advantage, that there are so many for free resources out there, that can help you to learn programming on your own. So there’s really no reason not to do so. And even if you just play around a bit, like I did, it will help you to face future challenges. Just always keep in mind, that it’s important to be comfortable with the basics in programming to solve harder problems. Everybody has to start small and if you do it for yourself, you don’t need to be perfect. Also nobody is perfect! 🙂

I hope you enjoyed this little story about my beginnings as a programmer. If you need some advise on this topic of starting to learn programming, you can just write me. So that I can make posts about topics, that might help you.
Until then, have a good one!

Please follow and like us:

The Fourier Spectrum Of An Image


Warning: count(): Parameter must be an array or an object that implements Countable in /homepages/10/d771928198/htdocs/clickandbuilds/RandomThoughtsScienceandProgramming/wp-content/plugins/kblog-metadata/kblog-author.php on line 332

Decomposing An Image Into Frequencies

The Fourier transform is a useful tool to disassemble a signal into its component frequencies. Amongst others it has applications in signal processing, quantum mechanics an even in Bioinformatics for the fast calculation of multiple alignments(10.1093/nar/gkf436).
So it’s no wonder, that it is also used in image processing. It is e.g. often used for the design of filters and for image convolutions in general, as it makes it possible to calculate those pretty fast. It can furthermore be used for the elimination of regular noise in images. But those are both topics for another time.


Today I will show you how to calculate the Fourier spectrum of an image, which is in any case useful to visualize what the Fourier transformation actually does. And you can also use it to investigate repeating features in your image. And from there it’s just one step further for using it for denoising an image.

How To Calculate It

I won’t got much into the theoretical background here right now, but I will give you the means to calculate it. If you want to get more intuition how the formulas for the Fourier transform are derived, I recommend you watching the following video of 3Brown1Blue.
To understand this part you should also at least know a little bit about the concept of imaginary numbers. But maybe I’ll do a quick tutorial in them at some point later. Just let me know, if you’d like to see that.

So what we actually need is the two dimensional discrete Fourier transformation, because pictures of course are two dimensional. Instead of discrete Fourier transformations there are also continuous ones, which can be applied to continuous functions.
So let me just give you the equations:

\(\hat{f}(p,q) = \frac{1}{\sqrt{n \cdot m}} \cdot \sum_{m = 0}^{M – 1}\sum_{n = 0}^{N – 1}f(n,m) \cdot e^{-\frac{i \cdot 2 \cdot \pi \cdot p \cdot m}{M}} \cdot e^{-\frac{i \cdot 2 \cdot \pi \cdot q \cdot n}{N}}\)

In this case N is the number of rows in the image and M the number of columns. The image itself is represented as the two dimensional function \(f(n, m)\) and \(\hat{f}(p,q)\) is the Fourier transformed image, where p and q are the new variables for indicating row and columns. And \(i\) is the imaginary unit. It’s a solution to the equation \(x^2 = -1\), which cannot be solved by any real number.

Now the value of each pixel is a complex number. Do you have any idea how we could visualize this in a sensible manner? Well… Complex numbers are basically just coordinates in a two dimensional space. So we could just take their absolute distance to the origin point in their coordinate system. Applying this to every pixel gives us the so-called Fourier spectrum of an image.

In the next parts I will step by step show you how to implement all of this and do some improvements on it. Luckily for us we don’t need to implement the Fourier transform ourselves. In R it is already contained in the form of the fft function and many other programming languages provide you with high performing implementations of it. Even though implementing the Fourier transform itself isn’t that big of a deal, performance could be. But maybe I’ll do some home-grown implementation of it at a later point in time.

Some Helper Functions

First we’ll need some helper functions. I was thinking about giving some of them their own blog post. But then again they’re not that interesting, so they would warrant that. We will also use some helper functions developed in my last post.

To make our life easier, a min and max function for images for be a nice, returning us the min and max grey-values of an image. The implementation is straight forward:

min.imageOneChannel <- function(object, na.rm = FALSE){
    return(min(object@imageMatrix, na.rm = na.rm))
}
max.imageOneChannel <- function(object, na.rm = FALSE){
    return(max(object@imageMatrix, na.rm = na.rm))
}

Next I’ll introduce the first of many point operations. Point operations change the value of one pixel independent of neighboring pixels. As the Fourier spectrum of an image might have values below zero or above one at certain pixels, we will need to rescale it so that it is inside the boundaries again. We can to this by with the following function:

affineRescale <- function(img){
    a <- max(img) - min(img)
    b <- min(img) / a

    return(affineGreyscaleTransformation(img = img, 1.0/a, -b))
}

After applying this function the minimal value will be equal to zero and the maximal equal to one and all other values will be scaled between them.

The Fourier Spectrum

So let’s implement the calculation of the Fourier spectrum of an image. You’ll see, that it is pretty straight-forward…

calculateFourierSpectrum <- function(img){

    imgFTrans <- fft(img@imageMatrix)

    ## calculate the eiclidean norm of the complex coordinates
    fourierSpectrum <- new("imageOneChannel", 
                           imageMatrix = sqrt(Re(imgFTrans) ^ 2 + 
                                              Im(imgFTrans) ^ 2))

    return(fourierSpectrum)
}

Basically just two/three (depending on if you count the return) lines of code and the nice thing about the native fft function in R, is that it keeps the dimensions of the matrix.

Now let’s load our today’s image

library(Raspository)
fence <- imageOneChannelFromJpeg("Fence.jpg")
Picture of a fence in Camphausen, Saarland.

It’s a fence, which has a lot of repeating structures. So it should give us some interesting result.

fftImg <- affineRescale(calculateFourierSpectrum(fence))
writeJPEG.imageOneChannel(fftImg, "fft1.jpg")
Fourier speactrum image transformation after affine rescaling

Huh?! What’s that? Nothing? Well, you can see a bit in the corners, but not much.

The Fourier Spectrum Scaled With The Logarithm

To make the spectrum better visible, we need to rescale it in a different manner. You remember, that the equation for the Fourier transform contains the exponent function? So what would probably be sensible would be to take the logarithm to rescale the image. But remember… The logarithm for smaller or equal to zero is not defined. So it’s a good idea to add a one before scaling with the logarithm. The log1p function in R will already do this for us.

logarithmicDynamicCompression <- function(img, c = NULL){
    if(is.null(c)){
        c <- 1 / log1p(max(img@imageMatrix))
    }

    return(new("imageOneChannel", imageMatrix = c * log1p(img@imageMatrix)))
}

So let’s try it again.

fftImg <- logarithmicDynamicCompression(calculateFourierSpectrum(fence))
writeJPEG.imageOneChannel(fftImg, "fft2.jpg")
Fourier speactrum image after logarithmic rescaling

Now that’s the Fourier spectrum of the fence. But there’s something not so nice about it. The origin of the spectrum is not in the middle, but on the borders. This is a certain property of the discrete form of the Fourier transformation. In the continuous one the origin lies in the middle.

Shifting The Fourier Spectrum Of An Image

But lucky for us there’s a quick fix for this minor problem. We just need to multiply the image with an image of the same size first. This second image contains one and minus ones arranged in a checkerboard like pattern. I won’t show you here, why this is so… But the Fourier transformation has many neat properties you can use, so that I’ll probably do a post about them in the future too.

For creating the checkerboard sign matrix I use data.table magic. Also probably doing a post or more posts about this awesome data structure in the future, as it makes of the life of a data scientist much easier.

library(data.table)
calculateFourierSpectrum <- function(img){


    ## generate checkerboard sign imageMatrix to shift it in the fourier domain
    DT <- CJ(1:nrow(img@imageMatrix), 1:ncol(img@imageMatrix))
    DT[, sign := (-1)^(V1 + V2)]
    shiftMatrix <- new("imageOneChannel", imageMatrix = as.matrix(dcast(DT, V1 ~ V2, value.var = "sign")[,-1]))

    shiftedImg <- img * shiftMatrix

    imgFTrans <- fft(shiftedImg@imageMatrix)


    ## calculate the eiclidean norm
    fourierSpectrum <- new("imageOneChannel", 
                           imageMatrix = sqrt(Re(imgFTrans)^2 + Im(imgFTrans)^2))

    return(fourierSpectrum)
}

Finally it’s time to also test that!

fftImg <- logarithmicDynamicCompression(calculateFourierSpectrum(fence))
writeJPEG.imageOneChannel(fftImg, "fft3.jpg")
Fourier speactrum image after logarithmic rescaling and shifting

And voila! Now I’m satisfied. 🙂
By the way the Fourier transform flips the direction of patterns. So the big vertical stripe could be the horizon. But it could also be a boundary artifacts… Something I’ll show you in more detail later. And also how to fix it.
All the stripes in the others direction represent frequencies across the picture. And what’s nice… You can change features of the image here in the Fourier space and then just as easily back-transform it.

You’ll definitely hear more from Fourier transformations from me. I hope you enjoyed this post. I’d love to hear from you.

Availability Of The Code

You can access a maintained form of the code in my GitHub repository Raspository under R/imageOneChannel, R/imagePointOperations and R/imageFeatures.

Please follow and like us:

Bibliography

Implementation Of Arithmetic Image Operators


Warning: count(): Parameter must be an array or an object that implements Countable in /homepages/10/d771928198/htdocs/clickandbuilds/RandomThoughtsScienceandProgramming/wp-content/plugins/kblog-metadata/kblog-author.php on line 332

Or More Object Oriented Programming (OOP) in R

As I told you in one of my recent posts I implemented a lot of functions related to image processing and built a lot on my image processing framework in the R programming language. Changes I did include methods I learned in the lecture, e.g. convolutions, filters of different kinds and calculation of statistical features. But I also did some code refactoring and included some new helper methods.
A set of those helper methods will be the topic of this blog-post. I’m talking about arithmetic image operators. We will need those a lot in later exercises. So it’s useful to have them implemented in an easy-to-use way.

So the basic idea ist just to overload the different arithmetic operators to make them useable with images of the type imageOneChannel1.
It should be possible to use the functions either with two images or one image and one scalar.

The is. Function

As a little helper we will implement the is.imageOneChannel function. Those functtions with the prefix is are used in R to test if any object is of a given type. It doesn’t save us much time from using the inherits function from the R base, but it makes the process at least somewhat more comfortable.
To be honest what I’ve implemented is just a glorifyed wrapper around the inherits function, but that’s generally how it is done. Of course that leaves us also with the possibility to write a little documentation using roxygen2. So here we go:

#' Is It A One Channel Image
#' 
#' @description Checks if an object is of type \code{\link[Raspository]{imageOneChannel}}.
#' 
#' @export
#'
#' @param x any R object.
#'
#' @return \code{TRUE}, if the object is of type \code{\link[Raspository]{imageOneChannel}}.
#'  \code{FALSE} otherwise.
#' 
#'
#' @examples
is.imageOneChannel <-function(x){
    return(inherits(x, "imageOneChannel"))
}

So, that was easy. And it’s not going to get any harder this post actually.

Arithmetic Image Operators

Overloading an arithmetic operator in R is pretty simple as well. Take a look at the following:

`+.imageOneChannel` <- function(A, B)

That will be the function header of our plus operator. Let’s analyse it bit by bit.
First the function name is sorrounded by backticks. That’s the way, how you define so called infix operators2, as they contain special symbols, which function names normally can’t.


Then we have the plus sign and after the dot the class name imageOneChannel, which indicates that this is the plus operator for this specific class. Later you don’t call this function as +.imageOneChannel, but just as plus. And as far I noticed it this special operator is always called, as long as one of the two parameters is of type imageOneChannel. But actually I don’t know what would happen, if the other summand would be of any other class, that also has its own plus operator defined. Anyone an idea? But I’ll probably test that at some point of time.
And last like in any other function <- function(A, B) assigns the function to the function name. As a plus function this contains just two parameters.

Now let me show you my full implementation of this arithmetic image operator with documentation:

#' Image Addition
#' 
#' @description Addition of one image with another image or a scalar. At least 
#' one of the two parameters has to be an image. If both are images, the operation
#' is executed entrywise.
#' 
#' @export
#'
#' @param A An image of class \code{\link[Raspository]{imageOneChannel}} or 
#' a scalar.
#' @param B An image of class \code{\link[Raspository]{imageOneChannel}} or 
#' a scalar.
#'
#' @return The resulting image of class \code{\link[Raspository]{imageOneChannel}}.
#' 
#'
#' @examples
`+.imageOneChannel` <- function(A, B){
    if(is.imageOneChannel(B)){
        if(is.imageOneChannel(A)){
            return(new("imageOneChannel", imageMatrix = A@imageMatrix + B@imageMatrix))
        }else{
            return(new("imageOneChannel", imageMatrix = A + B@imageMatrix))
        }

    }else{
        return(new("imageOneChannel", imageMatrix = A@imageMatrix + B))
    }

}

Pretty easy, isn’t it? We just need to check which one of the two parameters is an image and then access its matrix slot.
I mean for now our imageOneChannel also isn’t much more than a wrapper around a standard matrix in R. But that’s ok, because we definitly want to have some different behaviour defined later on than in a standard matrix. And We also need some methods, that not necessarily make sense for matrices containing other things than pixels.
Some people might know, that I’m not the biggest proponent of OOP and I actually think it’s overused. But if you have some data type, that should express some specific to it behaviour, it’s sensible. Just don’t implement a convoluted inheritance hierarchy with lots abd lots of stuff, you’ll never need and people need at least a PhD to understand. Keep it simple, stupid!

The other arithmetic image operators are implemented in the same manner. I’ll spare you the documentations this time, as they’re essentially all the same, but you can find them in my repository:

`-.imageOneChannel` <- function(A, B){
    if(is.imageOneChannel(B)){
        if(is.imageOneChannel(A)){
            return(new("imageOneChannel", imageMatrix = A@imageMatrix - B@imageMatrix))
        }else{
            return(new("imageOneChannel", imageMatrix = A - B@imageMatrix))
        }
    }else{
        return(new("imageOneChannel", imageMatrix  = A@imageMatrix - B))
    }

}
`*.imageOneChannel` <- function(A, B){
    if(is.imageOneChannel(B)){
        if(is.imageOneChannel(A)){
            return(new("imageOneChannel", imageMatrix = A@imageMatrix * B@imageMatrix))
        }else{
            return(new("imageOneChannel", imageMatrix = A * B@imageMatrix))
        }
    }else{
        return(new("imageOneChannel", imageMatrix = A@imageMatrix * B))
    }

}
`/.imageOneChannel` <- function(A, B){
    if(is.imageOneChannel(B)){
        if(is.imageOneChannel(A)){
            return(new("imageOneChannel", imageMatrix = A@imageMatrix / B@imageMatrix))
        }else{
            return(new("imageOneChannel", imageMatrix = A / B@imageMatrix))
        }
    }else{
        return(new("imageOneChannel", imageMatrix = A@imageMatrix / B))
    }

}

Testing The Arithmetic Image Operators

library(Raspository)

We can now do lot’s and lots of things with those operators. So first let’s load our test image…

kadser <- imageOneChannelFromJpeg("kadser.jpg")

Without looking at the original image let’s take a look at the negative of the picture.

negativeKadser <- 1 - kadser
writeJPEG.imageOneChannel(negativeKadser, "negativeKadser.jpg")
A negative of cat chilling in a weird position.

Spoopy, isn’t it?

And now let’s add this picture together with another one. We’ll have to average them, so that we don’t exceed the maximal grey-value.

owl0r <- imageOneChannelFromJpeg("owl0r.jpg")
addedImages <- (kadser + owl0r) / 2
writeJPEG.imageOneChannel(addedImages, "addedImages.jpg")
A cat chilling together with an owl.

Well, this added image won’t win us any contests, but at least we can say now, that we can add images together. That will be valuable later. You’ll see.

That’salso it for today. See you soon!

Availability Of The Code

You can access a maintained form of the code in my github repository Raspository under R/imageOperations.R and R/imageOneChannel.R.

Please follow and like us:

Rcpp – Integrating C++ Into R


Warning: count(): Parameter must be an array or an object that implements Countable in /homepages/10/d771928198/htdocs/clickandbuilds/RandomThoughtsScienceandProgramming/wp-content/plugins/kblog-metadata/kblog-author.php on line 332

Using Rcpp

Following my last post I will now show you how to integrate C++ code into a R package using Rcpp. Nothing revolutionary to do, but maybe it will still help some of you. 🙂
I will demonstrate this on my GitHub repository and R package BiocStyle::Githubpkg("David-J-R/Raspository"). And I will progress on optimizing the Dithering method I implemented last time. For the following I suppose, that you already use roxygen2.

Motivation

So at first you should ask yourself why you wanna introduce C++ code to your package. One very valid reason might be, that the C++ code is already there and you just want to integrate it into R or maybe even just write a wrapper in R around it.
The main reason why I did this was to speed up my code. Some operations in R are already pretty fast. But if you want to write a function with a loop that does complicated things, that would be torture to implement with apply functions, if even possible, C++ is a good choice.

There might of course be even more reasons, but we will save them for another time…
Let’s now get to how to do it.

Preparation For Integrating Rcpp Into R Packages

Before we can seamlessly integrate C++ using into our package we need to do some preparations. First we have to import the Rcpp package into our namespace and link to its dynamic library.
For this cause we add the following to the DESCRIPTION file of the package:

Imports: Rdpack, BEclear, futile.logger, stats, dequer, jpeg,    data.table, RobustRankAggreg, igraph, graphics, grDevices, stats, methods, abind, Rcpp
LinkingTo: Rcpp

And we also have to use the dynamic that will be created in our package. For this cause we add the following statement to the Raspository.R:

#' Raspository-package
#' 
#' @author David Rasp
#' 
#' @docType package
#' 
#' @import Rdpack
#' 
#' @title The Raspository
#' 
#' @useDynLib Raspository
"_PACKAGE"

As far as I understand it you could also add this statement anywhere into the documentation of your package.

Creating our function with Rcpp

So the translation of the dissipatePixel function from last time is pretty straight forward. It’s just changing the syntax a bit. I also added a roxygen skeleton. It also works the same as in R, just that you use two backslashes instead of a hash-tag to start a line. The namespace of Rcpp provides you with classes like NumericMatrix, that can be used similar to their R counterparts and the objects in R can be treated as them. If you want some more references, I would advise you to read the Rcpp quickref by Dirk Eddelbuettel. I have it open most of the time I’m working with it.

#include <Rcpp.h>
using namespace Rcpp;


//' Calculating the Dissipation of Pixels
//' 
//' @export dissipatePixel
//' @import Rcpp
//' 
//' 
//' @references \insertRef{Hagenburg2009}{Raspository}
//' 
//' @return the image as matrix after the next time step
//'
// [[Rcpp::export]]
SEXP dissipatePixel(const SEXP& imgOriginal, const double& minimalTreshold) {
    NumericMatrix img(imgOriginal);
    NumericMatrix imgAtNewStep(img.nrow(), img.ncol());

    for(std::size_t i = 0; i < img.nrow(); i++){
        for(std::size_t j = 0; j < img.ncol(); j++){

            double currentPixel = img(i,j);
            double dissipated = 0.0;

            if(currentPixel > 1.0 || currentPixel < minimalTreshold){

                double toDissipate = currentPixel; 

                if(currentPixel > 1){
                    toDissipate -= 1.0;
                }

                double tmp1 = toDissipate / 9.0;
                double tmp2 = toDissipate / 36.0;

                if(i > 0){
                    imgAtNewStep(i - 1, j) += tmp1;
                    dissipated += tmp1;
                }

                if(j > 0){
                    imgAtNewStep(i,j - 1) += tmp1;
                    dissipated += tmp1;
                }

                if(i < img.nrow() - 1){
                    imgAtNewStep(i + 1,j) += tmp1;
                    dissipated += tmp1;
                }

                if(j < img.ncol() - 1){
                    imgAtNewStep(i,j + 1) += tmp1;
                    dissipated += tmp1;
                }

                if( i > 0 && j > 0){
                    imgAtNewStep(i - 1,j - 1) += tmp2;
                    dissipated += tmp2;
                }

                if( i > 0 && j < img.ncol() - 1){
                    imgAtNewStep(i - 1,j + 1) += tmp2;
                    dissipated += tmp2;
                }

                if( i < img.nrow() - 1 && j > 0){
                    imgAtNewStep(i + 1,j - 1) += tmp2;
                    dissipated += tmp2;
                }

                if( i < img.nrow() - 1 && j > img.ncol() - 1){
                    imgAtNewStep(i + 1,j + 1) += tmp2;
                    dissipated += tmp2;
                }


            }else{

                double tmp1 = currentPixel / 9.0;
                double tmp2 = currentPixel / 36.0;

                if( i > 1 && img(i - 1,j) > currentPixel && img(i - 1,j) < 1){
                    imgAtNewStep(i - 1,j) += tmp1;
                    dissipated += tmp1;
                }

                if( j > 0 && img(i,j - 1) > currentPixel && img(i,j - 1) < 1){
                    imgAtNewStep(i,j - 1) += tmp1;
                    dissipated += tmp1;
                }

                if(i < img.nrow() - 1 && img(i + 1,j) > currentPixel && img(i + 1,j) < 1){
                    imgAtNewStep(i + 1,j) += tmp1;
                    dissipated += tmp1;
                }

                if(j < img.ncol() - 1 && img(i,j + 1) > currentPixel && img(i,j + 1) < 1){
                    imgAtNewStep(i,j + 1) += tmp1;
                    dissipated += tmp1;
                }

                if( i > 0 && j > 0 && img(i - 1,j - 1) > currentPixel && img(i - 1,j - 1) < 1){
                    imgAtNewStep(i - 1,j - 1) += tmp2;
                    dissipated += tmp2;
                }

                if( i > 0 && j < img.ncol() - 1 && img(i - 1,j + 1) > currentPixel && img(i - 1,j + 1) < 1){
                    imgAtNewStep(i - 1,j + 1) += 1.0 / 36.0 * currentPixel;
                    dissipated += 1.0 / 36.0 * currentPixel;
                }

                if( i < img.nrow() - 1 && j > 0 && img(i + 1,j - 1) > currentPixel && img(i + 1,j - 1) < 1){
                    imgAtNewStep(i + 1,j - 1) += tmp2;
                    dissipated += tmp2;
                }

                if( i < img.nrow() - 1 && j > img.ncol() - 1 && img(i + 1,j + 1) > currentPixel && img(i + 1,j + 1) < 1){
                    imgAtNewStep(i + 1,j + 1) += tmp2;
                    dissipated += tmp2;
                }

            }

            // add the non dissipated amount to the same pixel in next time-step
            imgAtNewStep(i,j) += currentPixel - dissipated;
        }
    }

return imgAtNewStep;
}

After writing and saving this function, we have to generate the RccpExports.cpp RccpExports.R files, which basically contains a wrapper for our C++ function, which then after compiling the code can be loaded in R.
We can comfortably do this with:

Rcpp::compileAttributes()

Then we should roxygenation our package, so that the namespace and the documentation get updated. We do this with:

roxygen2::roxygenise()

Now finally we can recompile the whole package and load the new version. The simplest way imo is to just use Ctrl+Shift+B in RStudio.

Performance comparison

Of we should also check, how much faster the new version is. I’ll use microbenchmark for this cause.

library(microbenchmark)
data <- matrix(runif(625), 25, 25)
microbenchmark(dissipatePixel(data, minimalTreshold = 0.05), dissipatePixelOld(data, minimalTreshold = 0.05))
## Unit: microseconds
##                                             expr      min       lq
##     dissipatePixel(data, minimalTreshold = 0.05)   70.891   72.368
##  dissipatePixelOld(data, minimalTreshold = 0.05) 4479.941 4622.600
##        mean   median        uq        max neval
##    95.23017   77.003   82.4445   1486.321   100
##  6514.79719 4677.444 4779.1125 161310.562   100

I think the result speaks for itself. This optimization was very worth while. The version implemented in C++ is about 60 times faster. That’s a lot. This could be the difference between one second and a minute. And of course longer calculations also suck more electricity, which means they have a larger CO2 fingerprint. So be nice to our atmosphere and code effective. 😉

That’s it for now.
Later,
David.

Availability Of The Code

You can access a maintained version of the code for the color spaces in my GitHub repository Raspository under src/dissipatePixel.cpp.

Please follow and like us:

Lattice Boltzmann Dithering


Warning: count(): Parameter must be an array or an object that implements Countable in /homepages/10/d771928198/htdocs/clickandbuilds/RandomThoughtsScienceandProgramming/wp-content/plugins/kblog-metadata/kblog-author.php on line 332

Revisiting Dithering

It’s been a while since my last post and an even longer while since the last one about image processing. However I’ll just catch up, where I left off. The last post in this category was about image dithering. A common usage of dithering is the conversion of a grey scale picture to a picture consisting just out of black and white pixels, saving space and also useful for some print media. On a further note dithered images have a nice visual style, I would call something like dollar bill style.

Conventional Dithering methods however also introduce artifacts to pictures and of course they also lose sharpness. So there’s the question, if one could approve on those conventional dithering methods in order to reduce artifacts and increase sharpness.
If you’re mathematically minded, you might also be worried, that conventional dithering is not rotationally invariant, which I assume is due to the fact that the dithering kernel is not symmetric and the way one iterates over the picture. Rotationally invariant basically means, that the output of a method doesn’t change, when you rotate your input image, which is of course a desirable trait.

So what about taking some inspiration from physics? After all physical laws of motion are rotationally invariant. Lattice Boltzmann methods are a way of computing dynamics in fluids. And they have been successfully employed for dithering by Hagenburg et al (10.1007/978-3-642-10520-3_91). The paper regarding this research also served as my main source for the method and its implementation, even though my specific implementation might differ in some details.
I’ll call the method Lattice Boltzmann Dithering during the rest of the article, as I didn’t find any other handy name for it.

The Idea Behind Lattice Boltzmann Dithering

The idea behind the method is to model the picture like a fluid with particles in it. In this case the grey values are the particles. A pixel that is black (with a grey value of 0) has so to say no particles in it. Those particles dissipate in the fluid over time. The time is modelled as discrete steps. They dissipate according to the following laws during each time step:

  • If a pixel has a value v greater than 1, where 1 is the maximal grey value, (v – 1) values are distributed amongst the neighbors of the pixel and the value of original pixel becomes 1.
  • If a pixel’s value v is lower than a minimal treshold, the user can set, v is distributed over its neighbors and the pixel itself becomes 0.
  • Otherwise a fraction of the pixels value is distributed to each neighbor, that has a larger value, while being smaller than 1. This fraction is subtracted from the pixel itself.

I hope you see, how with this laws pixels with a high value (many particles) attract more values. After some time steps you should more or less only be left with black and white pixels this way.

And the whole process is stopped, when the difference from one time step to the next one becomes sufficiently small. And that’s about it… Not really that complicated, if you ask me! 🙂
Of course the paper goes more into detail about the theoretical background and proves some things about the method. If you’re interested, you should absolutely read it.

Implementation of Lattice Boltzmann Dithering

And now let’s come to the part you all waited for… The implementation. First I’ll implement a method for the dissipation of the values from each pixel, which will e executed in each time step to get the next distribution of values.

One thing to keep in mind, when working with images is, that they have borders. And as you don’t want to dissipate the values across borders or worse access a part of the image, that isn’t there. So you have to treat the pixels at the borders differently. In the case of this method, I always check first, if a neighboring pixel is there or not.

dissipatePixel <- function(img, minimalTreshold = 0.01){


    imgAtNewStep <- matrix(c(0), ncol = ncol(img), nrow = nrow(img))

    for(i in seq(nrow(img))){
        for(j in seq(ncol(img))){
            if(img[i,j] > 1 || img[i,j] < minimalTreshold){

                toDissipate <- img[i,j] 

                if(img[i,j] > 1){
                    toDissipate <- toDissipate - 1
                }

                dissipated <- 0

                if( i > 1){
                    imgAtNewStep[i - 1,j] <- imgAtNewStep[i - 1,j] + 1.0 / 9.0 * toDissipate
                    dissipated <- dissipated + 1.0 / 9.0 * toDissipate
                }

                if( j > 1){
                    imgAtNewStep[i,j - 1] <- imgAtNewStep[i,j - 1] + 1.0 / 9.0 * toDissipate
                    dissipated <- dissipated + 1.0 / 9.0 * toDissipate
                }

                if(i < nrow(img)){
                    imgAtNewStep[i + 1,j] <- imgAtNewStep[i + 1,j] + 1.0 / 9.0 * toDissipate
                    dissipated <- dissipated + 1.0 / 9.0 * toDissipate
                }

                if(j < ncol(img)){
                    imgAtNewStep[i,j + 1] <- imgAtNewStep[i,j + 1] + 1.0 / 9.0 * toDissipate
                    dissipated <- dissipated + 1.0 / 9.0 * toDissipate
                }

                if( i > 1 && j > 1){
                    imgAtNewStep[i - 1,j - 1] <- imgAtNewStep[i - 1,j - 1] + 1.0 / 36.0 * toDissipate
                    dissipated <- dissipated + 1.0 / 36.0 * toDissipate
                }

                if( i > 1 && j < ncol(img)){
                    imgAtNewStep[i - 1,j + 1] <- imgAtNewStep[i - 1,j + 1] + 1.0 / 36.0 * toDissipate
                    dissipated <- dissipated + 1.0 / 36.0 * toDissipate
                }

                if( i < nrow(img) && j > 1){
                    imgAtNewStep[i + 1,j - 1] <- imgAtNewStep[i + 1,j - 1] + 1.0 / 36.0 * toDissipate
                    dissipated <- dissipated + 1.0 / 36.0 * toDissipate
                }

                if( i < nrow(img) && j > ncol(img)){
                    imgAtNewStep[i + 1,j + 1] <- imgAtNewStep[i + 1,j + 1] + 1.0 / 36.0 * toDissipate
                    dissipated <- dissipated + 1.0 / 36.0 * toDissipate
                }

                ## add the non dissipated amount to the same pixel in next time-step
                imgAtNewStep[i,j] <- imgAtNewStep[i,j] + (img[i,j] - dissipated)
            }else{

                dissipated <- 0
                currentPixel <- img[i,j]

                if( i > 1 && img[i - 1,j] > img[i,j] && img[i - 1,j] < 1){
                    imgAtNewStep[i - 1,j] <- imgAtNewStep[i - 1,j] + 1.0 / 9.0 * currentPixel
                    dissipated <- dissipated + 1.0 / 9.0 * currentPixel
                }

                if( j > 1 && img[i,j - 1] > img[i,j] && img[i,j - 1] < 1){
                    imgAtNewStep[i,j - 1] <- imgAtNewStep[i,j - 1] + 1.0 / 9.0 * currentPixel
                    dissipated <- dissipated + 1.0 / 9.0 * currentPixel
                }

                if(i < nrow(img) && img[i + 1,j] > img[i,j] && img[i + 1,j] < 1){
                    imgAtNewStep[i + 1,j] <- imgAtNewStep[i + 1,j] + 1.0 / 9.0 * currentPixel
                    dissipated <- dissipated + 1.0 / 9.0 * currentPixel
                }

                if(j < ncol(img) && img[i,j + 1] > img[i,j] && img[i,j + 1] < 1){
                    imgAtNewStep[i,j + 1] <- imgAtNewStep[i,j + 1] + 1.0 / 9.0 * currentPixel
                    dissipated <- dissipated + 1.0 / 9.0 * currentPixel
                }

                if( i > 1 && j > 1 && img[i - 1,j - 1] > img[i,j] && img[i - 1,j - 1] < 1){
                    imgAtNewStep[i - 1,j - 1] <- imgAtNewStep[i - 1,j - 1] + 1.0 / 36.0 * currentPixel
                    dissipated <- dissipated + 1.0 / 36.0 * currentPixel
                }

                if( i > 1 && j < ncol(img) && img[i - 1,j + 1] > img[i,j] && img[i - 1,j + 1] < 1){
                    imgAtNewStep[i - 1,j + 1] <- imgAtNewStep[i - 1,j + 1] + 1.0 / 36.0 * currentPixel
                    dissipated <- dissipated + 1.0 / 36.0 * currentPixel
                }

                if( i < nrow(img) && j > 1 && img[i + 1,j - 1] > img[i,j] && img[i + 1,j - 1] < 1){
                    imgAtNewStep[i + 1,j - 1] <- imgAtNewStep[i + 1,j - 1] + 1.0 / 36.0 * currentPixel
                    dissipated <- dissipated + 1.0 / 36.0 * currentPixel
                }

                if( i < nrow(img) && j > ncol(img) && img[i + 1,j + 1] > img[i,j] && img[i + 1,j + 1] < 1){
                    imgAtNewStep[i + 1,j + 1] <- imgAtNewStep[i + 1,j + 1] + 1.0 / 36.0 * currentPixel
                    dissipated <- dissipated + 1.0 / 36.0 * currentPixel
                }

                ## add the non dissipated amount to the same pixel in next time-step
                imgAtNewStep[i,j] <- imgAtNewStep[i,j] + (img[i,j] - dissipated)
            }
        }
    }



    return(imgAtNewStep)

}

Done that! Now the easy part…
But as the implementation in R with loops is incredibly inefficient I’ll just run 50 time steps this time. I will however implement this method at some later point in C++, where loops aren’t inefficient. This will also serve as a good example on how to integrate C++ code in R.

lbDithering <- function(img, epsilon, minimalTreshold){

    i <- 0
    while(TRUE){

        imgAtNewStep <- dissipatePixel(img = img, minimalTreshold = minimalTreshold)

        #if(norm(imgAtNewStep - img, type = "2") < epsilon){
        if(i >= 50){
            return(imgAtNewStep)
        }else{
            img <- imgAtNewStep
        }
        i <- i +1
    }
}

Usage

Now let’s reap the fruits of our coding and test the method on a picture I’ve selected. 🙂

birb <- imageBWFromJpeg("birb.jpg")
birbDithered <- lbDithering(birb@current, epsilon = 20, minimalTreshold = 0.05)
writeJPEG(round(birbDithered), "birbDithered.jpg")

Let me show you the resulting Jpeg:

A grainy picture of two pigeons lying down, dithered with the Lattice Boltzmann Method.
Image dithered with Lattice Boltzmann Method

Isn’t that beautiful? As comparison the original:

The original picture of two pigeons lying down.
Original image

And a conventionally dithered version:

birbDithered2 <- errorDiffusiondDithering(birb, method = "mea")
writeJPEG(birbDithered2@current, "birbDithered2.jpg")
A grainy picture of two pigeons lying down, dithered with the mea dithering.
Image dithered with minimized average error

You can see more structure in the one created with Lattice Boltzmann Dithering, don’t you? And also you can better understand how light and shadows are distributed.
So overall a pretty nice algorithm I would say! Although I like the dollar bill look of the conventional one as well.

So that’s it for now! Until soon…

Yours, David

Availability Of The Code

You can access a maintained version of the code for the color spaces in my GitHub repository Raspository under R/imageConversion.R.

And well… You might’ve noticed, that I used some of my methods/classes a bit different than the last few times. Having some distance from coding this stuff I noticed I have to change some things about it. Although I don’t know yet, if I’ll make a blog post about it or not. Do you wanna read one?

Please follow and like us:

Bibliography

D’Hondt Method


Warning: count(): Parameter must be an array or an object that implements Countable in /homepages/10/d771928198/htdocs/clickandbuilds/RandomThoughtsScienceandProgramming/wp-content/plugins/kblog-metadata/kblog-author.php on line 332

And Also Adams

In my last post I showed you the seat allocation method, called Sainte-Laguë/Schepers method. I recommend you reading it, before you continue with this post, as this post about the D’Hondt method heavily builds upon it.

Yepp, this vulture is cute!

There are actually two seat allocation methods, that are pretty similar to Sainte-Laguë/Schepers. One of them being the D’Hondt 1 method. The other one being the Adam’s method2.
Both are really remarkably similar to the method from my last post. With really the only difference being the way they round the seats. While Sainte-Laguë/Schepers uses standard round, D’Hondt uses the floor function. Meaning that it always round to the next lower integer. And Adams uses the ceiling function, which round to the next higher integer.

Now immediately you should scream: “STOP! What?! Adams uses the ceiling function?! So does this mean, that every party, that gets votes, gets at least one seat?!” The answer would be yesish. Yes, if your election doesn’t have an election threshold, every party, that gets votes, would at least get one seat.
“Well isn’t that incredibly unfavorable?” Yepp… But there are cases, where an allocation method like this could make sense. Not for regular election in my opinion, but for elections in parliaments. Let’s say, you have 20 mandates in a parliament, you want to distribute in a parliament with 300 elected politicians. Then the consideration, that it would be fair, if every party will get at least one mandate, could be made.

However the Adams method is incredibly uncommon. This Wikipedia article, which also served as my source for the method, only mentions the French parliament as example.

The D’Hondt method on the other side is pretty common. It’s actually the most common one in this year’s EU election.
And my source was also the corresponding German Wikipedia article.

Implementation D’Hondt Method

Luckily I don’t have to do much to implement those two methods. I just have to change a little bit about my function from last time.

seatAllocation <- function(votes, seats, roundMethod = round){
    ## calculate the initial divisor
    divisor <- roundMethod(sum(votes) / seats)

    ## get the initial seats per party
    seatsPerParty <- roundMethod(votes / divisor)

    ## if they already satisfy the seats to be assigned, return the seat allocation
    if(sum(seatsPerParty) == seats){
        return(seatsPerParty)
    }

    ## otherwise increment or decrement the divisor until 
    ## the result fits and then return it
    if(sum(seatsPerParty) < seats){
        while(sum(seatsPerParty) < seats){
            divisor = divisor - 1
            seatsPerParty <- roundMethod(votes / divisor)
        }
        return(seatsPerParty)
    }else{
         while(sum(seatsPerParty) > seats){
            divisor = divisor + 1
            seatsPerParty <- roundMethod(votes / divisor)
        }
        return(seatsPerParty = seatsPerParty)
    }

}

You see, what I did there? And why I love functional programming? Now by default, it’s the Sainte-Laguë/Schepers method and through giving the parameter roundMethod either the floor or ceiling function, we can make the D’Hondt and respectively Adams method out of it.
And we could even come up with some other rounding function in the future and use it.

The classic!

Test and Compare The Methods

And without further a due let’s test and compare the methods on our previous example.

votes <- c(AP = 11345, CVP = 563342, EP = 618713, OSP = 305952, PDP = 95001)
seatsSLS <- seatAllocation(votes, seats = 310, roundMethod = round)
seatsDH <- seatAllocation(votes, seats = 310, roundMethod = floor)
seatsA <- seatAllocation(votes, seats = 310, roundMethod = ceiling)

library(data.table)
DT <- rbind(data.table(party = names(seatsA), seats = seatsA, method = "Adams"), 
            data.table(party = names(seatsSLS), seats = seatsSLS, method = "Sainte-Laguë/Schepers"),
            data.table(party = names(seatsDH), seats = seatsDH, method = "D'Hondt"))

library(ggplot2)
g <- ggplot(DT, aes(x = party, y = seats, fill = method)) 
g <- g + geom_bar(stat = "identity", position = "dodge") 
g <- g + geom_text(aes(label=seats), position=position_dodge(width=0.9), vjust=-0.25)
g
Comparison of the different methods

Thanks, stackoverflow!
And you see… The actual difference isn’t big at all. The only thing one could say, is that Adams give a bonus to the small parties. D’Hondt method favors the big ones a bit. And Sainte-Laguë/Schepers is somehow in the middle.

And for me at least it’s really hard to say, which one is favorable. Sainte-Laguë/Schepers seems like a good compromise. However the differences more or less only affect small parties. But for them the difference is important. This doesn’t mean, that there’s no difference for large parties. On seat could mean the difference between majority and well… Not majority. Especially if you factor coalitions into the mix.
Maybe we will talk about possible problems in one of my next posts. I’m beginning to like this topic. I’m already thinking about becoming a lobbyist… lol.

Please follow and like us:

Sainte-Laguë/Schepers method


Warning: count(): Parameter must be an array or an object that implements Countable in /homepages/10/d771928198/htdocs/clickandbuilds/RandomThoughtsScienceandProgramming/wp-content/plugins/kblog-metadata/kblog-author.php on line 332

For Allocation of Seats in the EU Parliament

On Monday I had a talk over Discord with Boris Biba, who himself runs a blog. We wanted to do a cooperation for some time. The focus of his blog are philosophy and politics. And as I told him, that I’m interested in crunching numbers, the comming EU elections are the perfect opportunity for a cooperation.
First we talked about doing something regarding the Wahl-O-Mat. Now in hindsight it was probably good that we decided for something else, as the Wahl-O-Mat was taken offline just today.
Then Boris brought up that he wanted to a post about the seat allocation method, which is called Sainte-Laguë/Schepers method, for German votes in the EU election. And I thought to myself, that this is wonderful, as voting is basically a paradigm for statistics. So I would be able to implement a small algorithm.

So be also sure to check out the post, which you can find here, from Boris, if you’re able to read German!

More otter selfies means more creative common selfies! Like this one.

What I’ll be doing in this post, is explain you the seat allocation method called Sainte-Laguë/Schepers and then give you a demonstrative example for it. And as an easteregg I throw in some election posters for the imaginary parties, I’ll use in the example. I created those posters with Adobe Spark.

As a main source for my post, I took the corresponding article from the German Wahl-Lexikon.

Description of the Method

So there are basically three variants of this method, which all deliver the same result.
Two of them work by ranking the voting result. The other one by simple division, which is the one used for the German part of the EU election. It is either called iterative or divisor method.

The simple idea behind this divisor method is to find a divisor for the voting result, which delivers you the right amount of total seats, if you divide the voting results by it and then round them by standard rounding.

To find the right divisor, first the total amount of votes is divided by the number of seats to be assigned.

\(
divisor = \frac{\#votesTotal}{\#seats}
\)

The for each party the number of votes is divided by this divisor.

\(
seatsParty_{i} = \frac{\#votesParty_{i}}{divisor}
\)

And if the sum of the seats of all parties matches up with the amount to be assigned, we’re already done!
If not, we have to either increment or decrement the divisor depending on, if we have to few or to many seats.

Just think about that… If you increase the divisor, the amount of seats shrinks. And vice versa if you decrease the divisor, the amount of seats increases.

And so the divisor is adjusted and the final seats per party are obtained.

Who wouldn’t be proud of such a beautiful plumage?

Implementation of the Sainte-Laguë/Schepers method

And of course it wouldn’t be me, if I wouldn’t also implement the method.
Here we go…

seatAllocation <- function(votes, seats){
    ## calculate the initial divisor
    divisor <- round(sum(votes) / seats)

    ## get the initial seats per party
    seatsPerParty <- round(votes / divisor)

    ## if they already satisfy the seats to be assigned, return the seat allocation
    if(sum(seatsPerParty) == seats){
        return(list(divisor = divisor, seatsPerParty = seatsPerParty))
    }

    ## otherwise increment or decrement the divisor until 
    ## the result fits and then return it
    if(sum(seatsPerParty) < seats){
        while(sum(seatsPerParty) < seats){
            divisor = divisor - 1
            seatsPerParty <- round(votes / divisor)
        }
        return(list(divisor = divisor, seatsPerParty = seatsPerParty))
    }else{
         while(sum(seatsPerParty) > seats){
            divisor = divisor + 1
            seatsPerParty <- round(votes / divisor)
        }
        return(list(divisor = divisor, seatsPerParty = seatsPerParty))
    }

}

The function is basically the same as what I described under the last point in plain text. As always, if you have some questions or remarks regarding my implementation feel free to write me a comment!

Poor agenda or marketing gag? Hmmmm….

Example with the Sainte-Laguë/Schepers method

Now to test the method, let’s just come up with some arbitrary voting result for our imaginary parties introduced earlier. And of course plot them as a pie chart!

votes <- c(AP = 11345, CVP = 563342, EP = 618713, OSP = 305952, PDP = 95001)
votesSorted <- sort(votes, decreasing = TRUE)
names(votesSorted) <- paste0(names(votesSorted), " (", 
                               format(votesSorted/sum(votesSorted) * 100, digits = 2), "%)")
pie(votesSorted)
Pie chart of the voting result

Subsequently, let’s test what result the method delivers and if the percentages match up approximately.

result <- seatAllocation(votes, 310)

OK, first let’s visualize the result. But let’s not use a pie chart again. Because to be honest they can be misleading. This time we will use a waffle chart, which displays the actual seats.
Of course we also need to do some preprocessing. We want the parties ordered after their size and we won’t their percentage of seats in the legend.

seatsPerParty <- result$seatsPerParty
seatsPerParty <- sort(seatsPerParty, decreasing = TRUE)
names(seatsPerParty) <- paste0(names(seatsPerParty), " (", 
                               format(seatsPerParty/sum(seatsPerParty) * 100, digits = 2), "%)")
waffle::waffle(seatsPerParty)
Waffle diagram for the seat allocation

Well, there’s some difference in the percentage, but that’s to be expected as you can’t distribute fractions of seats between the parties.

Outlook

Of course there are many other methods for allocating seats in an election. Some that are equivalent to this one and others that are not. And if you’re interesting in them, I would encourage you to write me.
If you like, we can look at a bunch of them an then compare them. And we could also take a look at things like overhang seat or different kinds of voting. I think it’s a nice topic for making plots.

By the way if you also wanna read this post in German, check the following link out!

Please follow and like us:

Image Dithering


Warning: count(): Parameter must be an array or an object that implements Countable in /homepages/10/d771928198/htdocs/clickandbuilds/RandomThoughtsScienceandProgramming/wp-content/plugins/kblog-metadata/kblog-author.php on line 332

I showed you a lot of elementary stuff regarding image processing lately, so now it’s time to do something nice. What I will show you is called image dithering. But first let me give you a motivation for it

Let’s say you have a grayscale picture and you want it to paint with really only black and white. Therefore I will demonstrate you, what happens, if you just “round” the pixel values to black and white. For this purpose I will use a lot of functions I implemented in the previous posts. At first we have to load the picture.

library(Raspository)
imgColor <- imageRGBFromJpeg("Wagner.jpg")
plot(imgColor)
plot of chunk unnamed-chunk-1
A picture from the Nibelungen Halle in Königswinter of a monument to Richard Wagner, who had a neckbeard, before the invention of the internet.

Next step is, that we have to convert it to grayscale. Remember the function I implemented in this post?

img <- imageBWFromRGB(imgColor, c(0.3, 0.59, 0.11))
plot(img)
The same picture as a gray scale image.

Our eyes don’t treat all the colors equally, that’s why I choose this proportions for the colors to the grayscale. I took the values from this post on tutorialspoint.

Now let’s just quick an dirty convert this to completely black and white.

imgTmp <- img
imgTmp@current <- round(imgTmp@current)
plot(imgTmp)
The same picture as a with gray scales either rounded to black or white.

Wow! We certainly lost basically all the information in the picture.

Implementing Image Dithering

Before I show you how to implement a version of dithering I will shortly explain you the idea behind it. Let me ask you one question… Is there really such a thing as gray? Or how exactly would you define gray?
Quite simply as a mixture of black and white. Now think about colors and the screen, you’re probably sitting in front of. You have basically only three colors (RGB!). And each of the pixels on your screen consists of three sub-pixels, one of each of those three colors.
You perceive them not as individual dots, because their too close together for your eyes to distinguish them. Now’s the question, if we could something similar in this black and white case? And of course: Yes, we can! And it is called image dithering, which is by the way also applicable to colored images.

The idea now is that you iterate over all of your pixels ans apply to each of them still some kind of round function. But then you propagate the difference of the original pixel and the rounded pixels to its neighbors, that are still to be processed.
But of course there are also different methods in doing so. I’ll show you to of them today.

Floyd Steinberg dithering

Let’s begin with the Floyd Steinberg Algorithm. I suggest you to read the corresponding Wikipedia article, as it is very straightforward.

And my implementation is also pretty straightforward.

floydSteinbergDithering <- function(img, transformPaletteFunction = round){
    pixel <- img@current

    n <- dim(pixel)[1]
    m <- dim(pixel)[2]

    for(y in seq(m)){
        for(x in seq(n)){

            oldPixel <- pixel[x,y]
            newPixel <- transformPaletteFunction(oldPixel)

            error <- oldPixel - newPixel

            pixel[x,y] <- newPixel


                if(x < n){
                    pixel[x + 1, y] <- pixel[x + 1, y] + error * 7/16
                }

                if(x > 1 && y < m){
                    pixel[x - 1, y + 1] <- pixel[x - 1, y + 1] + error * 3/16
                }

                if(y < m){
                    pixel[x, y + 1] <- pixel[x, y + 1] + error * 5/16
                }

                if(x < n && y < m){
                    pixel[x + 1, y + 1] <- pixel[x + 1, y + 1] + error * 1/16
                }

        }
    }

    ditheredImage <- new(class(img)[[1]], original = img@original,
                         current = pixel, operations = img@operations)

    return(cropPixels(ditheredImage))
}

For the future some kind of Kernel function would be nice to be able to apply different kernels to pictures. But now let’s test it.

imgFS <- floydSteinbergDithering(img)
plot(imgFS)
The same picture with the Floyd-Steinberg algorithm applied to it.

That’s awesome! It almost looks like we had different gray values in our pictures. And there are just some minor artifacts introduced by it, meaning some appearent structures, that aren’t actually present in the original.

Now let’s try another method, which has a larger kernel.

Minimized Average Error Dithering

This dithering method was introduced by Jarvis et al from the famous Bell Lab in 1976. So you see that this whole field is pretty old. And some of you might remember a time, where it was actually difficult to transmit data from one location to another. I still remember being a six year old child waiting minutes on the NASA homepage to load one picture of a stellar nebular. Today image compression is of course still important for things like dynamic homepages, especially if they are mobile friendly.

OK, now let’s come to the actual method. It is called minimized average error. And again the Wikipedia article on it is pretty good.

This time the neighborhood of your pixel is increased by a range of one. Let me show you the implementation…

minimizedAverageErrorDithering <- function(img, transformPaletteFunction = round){
    pixel <- img@current

    n <- dim(pixel)[1]
    m <- dim(pixel)[2]

    for(y in seq(m)){
        for(x in seq(n)){

            oldPixel <- pixel[x,y]
            newPixel <- transformPaletteFunction(oldPixel)

            error <- oldPixel - newPixel

            pixel[x,y] <- newPixel

                if(x < n){
                    pixel[x + 1, y    ] <- pixel[x + 1, y    ] + error * 7/48
                }
                if(x < n - 1){
                    pixel[x + 2, y    ] <- pixel[x + 2, y    ] + error * 5/48
                }

                if(x > 2 && y < m){
                    pixel[x - 2, y + 1] <- pixel[x - 2, y + 1] + error * 3/48
                }

                if(x > 1 && y < m){
                    pixel[x - 1, y + 1] <- pixel[x - 1, y + 1] + error * 5/48
                }

                if(y < m){
                    pixel[x    , y + 1] <- pixel[x    , y + 1] + error * 7/48
                }

                if(x < n && y < m){
                    pixel[x + 1, y + 1] <- pixel[x + 1, y + 1] + error * 5/48
                }

                if(x < n - 1 && y < m){
                    pixel[x + 2, y + 1] <- pixel[x + 2, y + 1] + error * 3/48
                }

               if(x > 2 && y < m - 1){
                    pixel[x - 2, y + 2] <- pixel[x - 2, y + 2] + error * 1/48
                }

                if(x > 1 && y < m - 1){
                    pixel[x - 1, y + 2] <- pixel[x - 1, y + 2] + error * 3/48
                }

                if(y < m - 1){
                    pixel[x    , y + 2] <- pixel[x    , y + 2] + error * 5/48
                }

                if(x < n && y < m - 1){
                    pixel[x + 1, y + 2] <- pixel[x + 1, y + 2] + error * 3/48
                }

                if(x < n - 1 && y < m - 1){
                    pixel[x + 2, y + 2] <- pixel[x + 2, y + 2] + error * 1/48
                }
        }
    }

    ditheredImage <- new(class(img)[[1]], original = img@original,
                         current = pixel, operations = img@operations)

    return(cropPixels(ditheredImage))
}

You wanna see it’s effect, don’t you? Here you go…

imgMea <- minimizedAverageErrorDithering(img)
plot(imgMea)
The same picture with Minimized Average Error applied to it.

Do you see the difference? I think we got rid of the artifacts! Isn’t that amazing? I really love how demonstrative image processing is.
But that’s it for today… See you soon!

Availability Of The Code

You can access a maintained version of the code for the color spaces in my GitHub repository Raspository under R/imageConversion.R.

Please follow and like us:

Converting HSV To RGB


Warning: count(): Parameter must be an array or an object that implements Countable in /homepages/10/d771928198/htdocs/clickandbuilds/RandomThoughtsScienceandProgramming/wp-content/plugins/kblog-metadata/kblog-author.php on line 332

This post is now the sequel to this week’s post about color spaces. I think it will be a rather short post, where I just implement the conversion back HSV to RGB and then test my methods for demonstration.
In doing so I used the corresponding wikipedia post as a rough guideline.

So let me show you how I did it.

library(Raspository)
calculateRGBvalue <- function(H, C, X, m){
   if(H >= 0 && H <= 1){
       return(c(m + C, m + X, m))
   }else if(H >= 0 && H <= 2){
       return(c(m + X, m + C, m))
   }else if(H >= 0 && H <= 3){
       return(c(m, m + C, m + X))
   }else if(H >= 0 && H <= 4){
       return(c(m, m + X, m + C))
   }else if(H >= 0 && H <= 5){
       return(c(m + X, m, m + C))
   }else if(H >= 0 && H <= 6){
       return(c(m + C, m, m + X))
   }else{
       return(c(0,0,0))
   }
}
require(abind)
## Loading required package: abind
hsvArrayToRgb <- function(hsvArray){

    # Calculate the chroma
    C <- hsvArray[,,3] * hsvArray[,,2]

    H<- hsvArray[,,1] / 60

    X <-  C * (1 - abs(H %% 2 - 1))

    m <- hsvArray[,,3] - C
    rgb<-mapply(FUN = calculateRGBvalue, H = H, C = C, X = X, m = m)

    rgbArray<-abind(matrix(rgb[1,], nrow = nrow(hsvArray)),
                matrix(rgb[2,], nrow = nrow(hsvArray)),
                matrix(rgb[3,], nrow = nrow(hsvArray)),
                along = 3)

    return(rgbArray)
}


imageRGBFromHSV <- function(img){
    return(new("imageRGB", original = hsvArrayToRgb(img@original),
               current = hsvArrayToRgb(img@current),
               operations = img@operations))
}

So far so good… That’s basically just an one to one implementation of the method from the Wikipedia article with the one difference being, that I add m beforehand and not later.

So let’s come to the fun part!

image <- imageRGBFromJpeg("Mountain.jpg")
## Loading required package: jpeg
plot(image)
The original pictures

What a ginormous mountain! But it isn’t even the tallest mountain in Saarland, believe it! Now we can also plot this as an black and white picture.

plot(imageBWFromRGB(image))
Picture in black and white

Then let’s convert to HSV and plot the S (saturation) and V (value/brightness) channel as black and white picture.

hsv <- imageHSVFromRGB(image)
plot(as.raster(hsv@current[,,2]))
The saturation component of the pictures

How that looks like something from the 90s!

plot(as.raster(hsv@current[,,3]))
The brightness component of the pictures

OK, and now compare this to the conversion to black and white from RGB. It’s a worse picture in my opinion, but in comparison to the other one it holds the accurate information about the lighting condition of each pixel. You see, how this could be useful in Computer Vision e.g.?
Let’s just say, you’re searching for something in your pictures, that absorbs a lot of light. Then this representation could be of use.

And last, but not least we have to convert the image back to confirm that we actually get back the original picture.

plot(imageRGBFromHSV(hsv))
The pictures first converted to HSV and then back to RGB

Looks the same! I guess that concludes this blog post. And I already have some plans for the next few posts.
So, hopefully see you! 🙂

Please follow and like us:

Introduction To Color Space


Warning: count(): Parameter must be an array or an object that implements Countable in /homepages/10/d771928198/htdocs/clickandbuilds/RandomThoughtsScienceandProgramming/wp-content/plugins/kblog-metadata/kblog-author.php on line 332

For now we have only looked at black and white pictures, but of course I know y’all are craving for me talking about color pictures! And that’s, what I will do today… After I visited a lecture about it.
What I will exactly do today is implement a S4 class for colored pictures with an RGB color space similar to one for black and white picture and then implement one with an HSV-color space and add corresponding transformation methods between the two classes.
So pretty basic stuff as well, but especially the HSV version could be pretty useful for manipulating pictures or further using them in Computer Vision procedures.

A colorful picture with the upper half being the sky with clouds and the lower half the view over the Spandau Citadel and waters in the background.
A picture taken in Berlin-Spandau at the Citadel – Scenes like this are a nice remainder would a beautiful thing color vision is. 🙂

RGB Color Space

Before we have only looked at pixels with a one dimensional value between zero and one.
In a more mathematical way, we looked at a function:

\(
\mathbb{N} \times \mathbb{N} \rightarrow [0,1]
\)

If this notation confuses you… Just ignore it! But in other words: The picture is a function of two natural numbered coordinates, that gives you back a real numbered value (an intensity) between zero and one.

If we want to do the same with colored images it is useful to define each pixel as having three dimension. One for Red, one for Green and one for Blue. Therefore the name RGB.
This way we would have a function:

\(
\mathbb{N} \times \mathbb{N} \rightarrow [0,1] \times [0,1] \times [0,1]
\)

Again in other words: The picture is a function of two natural numbered coordinates, that gives you back a tuple1 of three real numbered values between zero and one.
By the way this RGB color space is probably what you’re looking at right now, as computer monitors usually use it, but of course things in reality are more complicated. There are more than three wavelengths of light, but as our color vision system has only three different color receptors, a RGB color space is good enough for the human eye. But of course in scientific system, it could very well be, that you want to have more than three colors in a picture. Coming originally from Molecular Biology I know that there are many distinct wavelengths, that have some Biological importance. E.g. where certain molecules or functional groups of molecules absorb light.
So you could imagine using color spaces for higher dimension for some scientific applications.

But let’s return to the topic of RGB color spaces. In programming terms one way of implementing it would be with a three dimensional array, which is pretty straight forward.
So without further ado let me show you my implementation. It’s pretty simple, as the jpeg already returns the right kind of array, when reading a jpeg.

imageRGB <- setClass("imageRGB", slots=list(original="array", current="array", operations="list"))
imageRGBFromJpeg <-function(pathToJpeg){
    require(jpeg)
    image<- readJPEG(pathToJpeg)
    return(new("imageRGB", original = image, current = image, operations = list()))
}
plot.imageRGB <- function(object){
    plot(as.raster(object@current))}

Converting To Black And White

Now I also want a method to convert an RGB image into a black and white one.

library(Raspository)
imageBWFromRGB <- function(img, chPortion = c(0.33, 0.33, 0.33)){
    if(sum(chPortion) > 1){
        stop("Channel portions mustn't add up to more than one.")
    }
    original <- img@original[,,1] * chPortion[1] + img@original[,,2] * chPortion[2] + img@original[,,3] * chPortion[3]
    current <- img@current[,,1] * chPortion[1] + img@current[,,2] * chPortion[2] + img@current[,,3] * chPortion[3]
    return(new("imageBW", original = original, current = current, operations = img@operations))
}

The parameter chPortion is a vector that indicates how much each channel from the original color space should contribute to the result black and white picture. As a matter of course its components aren’t allowed to add up to more than one.

HSV Color Space

Another color space I want to show you today is the HSV. While you could imagine RGB is a cube, in the same analogy HSV would be a cylinder. With one polar coordinate H, which stands for hue. This coordinate determines the color at a given pixel.
The other two coordinates are S saturation and V the value (brightness) of a pixel. Those two coordinates are rational numbers between zero and one as before.

Now you could ask, why this representation could be useful and the answers are many. But let me just name you the editing of pictures as one of them. Let’s say you have a picture that is to dark, but the colors itself are OK. Then you could use this HSV color space for changing the value/brightness, while leaving hue and saturation the same.
But of course there are many other color spaces having other neat properties you might wanna exploit.
However instead of dreaming about the perfect color space, which feels like magic, let’s just implement this one, including the conversion from RGB.
We’ll also need some helper functions for it. And please, if you don’t understand something, just ask me in the comments. There should be no shame in don’t understanding something and I will be happy to give you an explanation to help you understanding it.

imageHSV <- setClass("imageHSV", slots=list(original="array", current="array", operations="list"))

calculateHUE <- function(max, maxIndex, min, r, g, b){
    h <- 0.0
    if(max == min){
        return(h)
    }else if(maxIndex == 1){
        h <- 60.0 * ((g - b)/(max - min)) 
    }else if(maxIndex == 2){
        h <- 60.0 * (2.0 + (b - r)/(max - min)) 
    }else if(maxIndex == 3){
        h <- 60.0 * (4.0 + (r - g)/(max - min)) 
    }

    # if the value is negativ add 360° to it
    if(h >= 0){
        return(h)
    }else{
        return(h + 360)
    }
}

rgbArrayToHsv <- function(rgbArray){
    # get the maximal color and its index in each pixel
    max <- apply(rgbArray, c(1,2), max)
    maxIndex <-apply(rgbArray, c(1,2), which.max)
    # get the minimal color in each pixel
    min <- apply(rgbArray, c(1,2), min)

    # calculate the hue for each pixel
    h <- mapply(FUN = calculateHUE, max = max, maxIndex = maxIndex, min = min, 
                r = rgbArray[,,1], g = rgbArray[,,2], b = rgbArray[,,3])
    # convert vector back to matrix
    h <- matrix(h, ncol = ncol(max))

    # calculate saturation
    s <- (max - min)/max
    # set values to zero, where max is 0 (division by zero -> NA)
    s[is.na(s)] <- 0
    # max is equal to v (value/brightness)
    v <- max

    # bind matrices together to array and return
    require(abind)
    hsvArray <- abind(h, s, v, along = 3)
    return(hsvArray)
}

imageHSVFromRGB <- function(img){
    return(new("imageHSV", original = rgbArrayToHsv(img@original), 
               current = rgbArrayToHsv(img@current), 
               operations = img@operations))
}

Phew, this took me more effort than expected. Luckily I feel comfortable using the various apply functions in R. They might be alien at first, but there’s a lot of neat stuff you can do, if you can use them.
But one side note… The calculateHUE function is very time consuming, because it is called a lot. So I think it would be worth while to re-implement it at some point of time in C++ to speed up things. Probably a good exercise for the future. I’ll keep you posted. 🙂

Anyway that’s enough for me today… David tired now. 😮 The conversion from HSV to RGB will probably come the next few days.
So, see ya!

Availability Of The Code

And of course you can access a maintained version of the code for the color spaces in my GitHub repository Raspository under R/imageRGB.R and R/imageHSV.R.

Please follow and like us: