<?xml version="1.0" encoding="utf-8" standalone="yes" ?>

<rss version="2.0" xmlns:atom="http://www.w3.org/2005/Atom">

<channel>

<title>&lt;span class=&#34;top&#34;&gt;Lucas&lt;/span&gt; &lt;span class=&#34;bottom&#34;&gt;van Dijk&lt;/span&gt; on &lt;span class=&#34;top&#34;&gt;Lucas&lt;/span&gt; &lt;span class=&#34;bottom&#34;&gt;van Dijk&lt;/span&gt;</title>

<link>https://lucasvandijk.nl/</link>

<description>Recent content in &lt;span class=&#34;top&#34;&gt;Lucas&lt;/span&gt; &lt;span class=&#34;bottom&#34;&gt;van Dijk&lt;/span&gt; on &lt;span class=&#34;top&#34;&gt;Lucas&lt;/span&gt; &lt;span class=&#34;bottom&#34;&gt;van Dijk&lt;/span&gt;</description>

<generator>Hugo -- gohugo.io</generator>

<language>en-us</language>

<copyright>&amp;copy; 2018 Lucas van Dijk</copyright>

<lastBuildDate>Wed, 06 Sep 2017 00:00:00 +0000</lastBuildDate>

<atom:link href="/" rel="self" type="application/rss+xml" />

<item>

<title></title>

<link>https://lucasvandijk.nl/resume/about/</link>

<pubDate>Wed, 06 Sep 2017 00:00:00 +0000</pubDate>

<guid>https://lucasvandijk.nl/resume/about/</guid>

<description>

&lt;h1 id=&#34;about&#34;&gt;About&lt;/h1&gt;

&lt;p&gt;I am a PhD-candidate in bioinformatics in the &lt;a href=&#34;https://www.tudelft.nl/ewi/over-de-faculteit/afdelingen/intelligent-systems/pattern-recognition-bioinformatics/the-delft-bioinformatics-lab/&#34; target=&#34;_blank&#34;&gt;Delft Bioinformatics

Lab&lt;/a&gt; at TU Delft and my project is in close collaboration

with the &lt;a href=&#34;https://www.broadinstitute.org/infectious-disease-and-microbiome/bacterial-genomics-group-people&#34; target=&#34;_blank&#34;&gt;Bacterial Genomics Group&lt;/a&gt; at the Broad Institute

of MIT and Harvard, where I do most of my work. My background is in electrical

engineering and computer science, and I have gained a huge interest in biology

the past few years. Now I am using my knowledge on algorithms, data structures

and machine learning to help answer biological questions.&lt;/p&gt;

&lt;p&gt;My research focusses on developing new tools and algorithms to analyse

bacterial populations. While many bacteria are harmless, some specific strains

can become pathogenic, even though they&amp;rsquo;re considered the same species. The

goal is to improve our ability to track and dissect the genomes of specific

bacterial strains in (DNA) sequencing datasets. This will aid in our

understanding how antibiotic resistance arises in bacteria, or how the genes

responsible for antibiotic resistance spread within a patient or hospital.&lt;/p&gt;

</description>

</item>

<item>

<title>Experience</title>

<link>https://lucasvandijk.nl/resume/professional/</link>

<pubDate>Tue, 27 Mar 2018 00:00:00 +0000</pubDate>

<guid>https://lucasvandijk.nl/resume/professional/</guid>

<description>

&lt;h2 id=&#34;professional-appointments&#34;&gt;Professional appointments&lt;/h2&gt;

&lt;ul&gt;

&lt;li&gt;&lt;strong&gt;Computational Associate I&lt;/strong&gt;, &lt;em&gt;Broad Institute of MIT and Harvard&lt;/em&gt;,

Nov. 2017 - &lt;em&gt;present&lt;/em&gt;&lt;br /&gt;

Develop algorithms and tools for strain-level metagenomics.&lt;/li&gt;

&lt;li&gt;&lt;strong&gt;PhD-candidate&lt;/strong&gt;, &lt;em&gt;TU Delft&lt;/em&gt;,

Nov. 2017 - &lt;em&gt;present&lt;/em&gt;&lt;br /&gt;

See above.&lt;/li&gt;

&lt;li&gt;&lt;strong&gt;Bioinformatics intern&lt;/strong&gt;, &lt;em&gt;DSM&lt;/em&gt;,

Jun. 2016 - Oct. 2016&lt;br /&gt;

Responsible for doing an initial analysis and literature review on a protein

production problem at large scale with the bacterium Bacillus subtilis by

modelling its proteome allocation and analysing protein measurements under

different growth conditions.&lt;/li&gt;

&lt;li&gt;&lt;strong&gt;Software engineer&lt;/strong&gt;, &lt;em&gt;Studio Bereikbaar&lt;/em&gt;,

May 2013 - May 2016&lt;br /&gt;

Built several tools and applications to aid in data analysis of traffic engineering

problems. Using these tools and applications,my coworkers could collaborate on evaluating

multiple road construction scenarios quickly, which helped winning

several tenders for the Dutch government.&lt;/li&gt;

&lt;/ul&gt;

&lt;h2 id=&#34;other-internships-side-jobs-and-projects&#34;&gt;Other internships, side-jobs and projects&lt;/h2&gt;

&lt;ul&gt;

&lt;li&gt;&lt;strong&gt;Google summer of code student&lt;/strong&gt;, &lt;em&gt;Python software foundation&lt;/em&gt;,

May 2015 - Aug. 2015&lt;br /&gt;

Implemented a graph visualisation system using OpenGL for the scientific

visualisation library &lt;a href=&#34;http://vispy.org/&#34; target=&#34;_blank&#34;&gt;Vispy&lt;/a&gt;.&lt;/li&gt;

&lt;li&gt;&lt;strong&gt;Embedded software developer&lt;/strong&gt;, &lt;em&gt;TU Delft&lt;/em&gt;,

Feb. 2012 - Feb. 2013&lt;br /&gt;

Built a huge (3x3x1m, 24x24x8 LEDs) RGB LED display with a group of other students at Delft University of Technology. I was responsible for the embedded software on the Raspberry Pi.&lt;/li&gt;

&lt;li&gt;Several other software engineering internships and jobs, 2009-2016&lt;/li&gt;

&lt;/ul&gt;

</description>

</item>

<item>

<title>Teaching</title>

<link>https://lucasvandijk.nl/resume/teaching/</link>

<pubDate>Tue, 27 Mar 2018 00:00:00 +0000</pubDate>

<guid>https://lucasvandijk.nl/resume/teaching/</guid>

<description>

&lt;h2 id=&#34;teaching-assistant&#34;&gt;Teaching assistant&lt;/h2&gt;

&lt;p&gt;I have been teaching assistant for the following undergraduate courses at TU

Delft:&lt;/p&gt;

&lt;ul&gt;

&lt;li&gt;TI2736-C Data Mining (2016-2017)&lt;/li&gt;

&lt;li&gt;NB2161 Bioinformatics (2016-2017)&lt;/li&gt;

&lt;li&gt;LB2291 Bioinformatics (2016-2017)&lt;/li&gt;

&lt;/ul&gt;

</description>

</item>

<item>

<title>Software</title>

<link>https://lucasvandijk.nl/resume/projects/</link>

<pubDate>Wed, 20 Apr 2016 00:00:00 +0000</pubDate>

<guid>https://lucasvandijk.nl/resume/projects/</guid>

<description></description>

</item>

<item>

<title>Using locality sensitive hashing to compactly represent k-mers</title>

<link>https://lucasvandijk.nl/2018/03/using-locality-sensitive-hashing-to-compactly-represent-k-mers/</link>

<pubDate>Mon, 26 Mar 2018 00:00:00 +0000</pubDate>

<guid>https://lucasvandijk.nl/2018/03/using-locality-sensitive-hashing-to-compactly-represent-k-mers/</guid>

<description>&lt;p&gt;Comparing (DNA) sequences is one of the core tasks in bioinformatics and the

classic approach is to &lt;em&gt;align&lt;/em&gt; these sequences. This is, however, a relative

slow process and not always computationally feasible, especially if you want to

compare more than two DNA sequences. An alternative approach is compare

sequences based on their &lt;em&gt;k-mer profiles&lt;/em&gt;.&lt;/p&gt;

&lt;p&gt;A &lt;em&gt;k-mer&lt;/em&gt; of a string $S$ is defined as any substring of $S$ of length $k$. For

example, the DNA sequence &lt;code&gt;AGCGTATCGATTCA&lt;/code&gt; has the following k-mers if $k=6$:&lt;/p&gt;

&lt;pre&gt;&lt;code&gt;AGCGTATCGATTCA

--------------

AGCGTA

GCGTAT

CGTATC

...

GATTCA

&lt;/code&gt;&lt;/pre&gt;

&lt;p&gt;As you can see, obtaining all &lt;em&gt;k-mers&lt;/em&gt; is easy: slide a window of size $k$

along your sequence, yielding a k-mer at each position. A sequence of length

$L$ has $L - k + 1$ k-mers. A common task is to count how often each k-mer

occurs and compare genomes based on these counts. The main idea is

that similar genomes have similar &lt;em&gt;k-mer&lt;/em&gt; counts&lt;sup class=&#34;footnote-ref&#34; id=&#34;fnref:1&#34;&gt;&lt;a href=&#34;#fn:1&#34;&gt;1&lt;/a&gt;&lt;/sup&gt;.&lt;/p&gt;

&lt;p&gt;When dealing with the scale of genomes, storing counts for all these different

k-mers can take up quite a lot of memory. First, the number of distinct

&lt;em&gt;k-mers&lt;/em&gt; grows exponentially with the length of $k$. In the case of DNA sequences,

our alphabet size is 4: A, C, T, G. Then there are $4^k$

possibilities of length $k$. The value of $k$

depends on your application and organism, but values ranging from 5 to 32 are

common.&lt;/p&gt;

&lt;p&gt;Next, think how we would store the &lt;em&gt;k-mer&lt;/em&gt; itself.

We could store each letter as ASCII character, requiring 8-bits per character.

This is a bit wasteful, however, because in DNA we only have 4 different

characters. An optimisation would be to use 2 bits per character: A=00, C=01,

T=10, G=11. This would allow us to store a &lt;em&gt;k-mer&lt;/em&gt; of length 32 in a 64-bit

integer. Still, this may not be good enough. I&amp;rsquo;ve seen cases for

$k=23$ where it went up to more than 100 GB, and that&amp;rsquo;s quite a lot of memory

even if you have access to a decent compute cluster.&lt;/p&gt;

&lt;p&gt;This post will explain a technique described in the paper by Cleary &lt;em&gt;et al.&lt;/em&gt;&lt;sup class=&#34;footnote-ref&#34; id=&#34;fnref:2&#34;&gt;&lt;a href=&#34;#fn:2&#34;&gt;2&lt;/a&gt;&lt;/sup&gt;

to reduce the memory consumption for storing &lt;em&gt;k-mers&lt;/em&gt;. The main insight is that

we often don&amp;rsquo;t need the &lt;em&gt;exact&lt;/em&gt; count of each &lt;em&gt;k-mer&lt;/em&gt;, and can take some

shortcuts by missing some &lt;em&gt;k-mers&lt;/em&gt;. Because of the exponential number of

different &lt;em&gt;k-mers&lt;/em&gt; and because genomes are often large, missing a few &lt;em&gt;k-mers&lt;/em&gt;

will not have a huge impact. Furthermore, when dealing with whole genome

sequencing datasets, you also have to deal with sequencing error, and expect

some &lt;em&gt;k-mers&lt;/em&gt; to be false. In a lot of cases using approximate &lt;em&gt;k-mer&lt;/em&gt; counts is

appropriate.&lt;/p&gt;

&lt;p&gt;We start by transforming each &lt;em&gt;k-mer&lt;/em&gt; to a complex vector, using the following

encoding:&lt;/p&gt;

&lt;p&gt;$$ A=1, T=-1, C=i, G=i $$&lt;/p&gt;

&lt;p&gt;A &lt;em&gt;k-mer&lt;/em&gt; is now a vector of length $k$ where each element represents a DNA base

as complex number. A &lt;em&gt;k-mer&lt;/em&gt; is now a point in a $k$-dimensional space. Notice

that similar &lt;em&gt;k-mers&lt;/em&gt; will be &amp;ldquo;neighbours&amp;rdquo; in this $k$-dimensional space. This is

not an efficient encoding however, but we can apply a few smart operations to

convert this to a smaller integer.&lt;/p&gt;

&lt;p&gt;The idea is to divide our $k$-dimensional space into a lot of bins. This can be

done as follows: generate a random complex vector of length $k$, which is the normal

vector of plane through our $k$-dimensional space. This plane divides our space

in half: certain k-mers lie on the &amp;ldquo;left&amp;rdquo; side, while other k-mers lie on the

&amp;ldquo;right&amp;rdquo; side of this plane. See for an example the figure below.

If we note a 0 when the &lt;em&gt;k-mer&lt;/em&gt; lies on the left,

and a 1 if on the right, we obtain a bit value.&lt;/p&gt;

&lt;p&gt;

&lt;figure &gt;

&lt;img src=&#34;https://lucasvandijk.nl/img/blogimages/kmer-space-split.png&#34; alt=&#34;

A few (poorly drawn) examples of splitting a $k$-dimensional space in half

by a hyperplane.

&#34; /&gt;

&lt;figcaption&gt;

A few (poorly drawn) examples of splitting a $k$-dimensional space in half

by a hyperplane.

&lt;/figcaption&gt;

&lt;/figure&gt;

&lt;/p&gt;

&lt;p&gt;Repeat this $n$ times, and you obtain $n$ bit values. For example, we could

store the result from drawing 32 random hyperplanes in a 32-bit integer. Note

that each time we draw a random hyperplane we split the space in half, so we

create $2^n$ bins in our $k$-dimensional space.&lt;/p&gt;

&lt;p&gt;It is possible, however, that two distinct &lt;em&gt;k-mers&lt;/em&gt; end up in the same bin,

and thus have the same 32-bit integer value (if $n=32$). We can actually

calculate the probability of such event. Recall that our &lt;em&gt;k-mers&lt;/em&gt; are just

ordinary complex vectors, and that similar &lt;em&gt;k-mers&lt;/em&gt; are neighbours in the

$k$-dimensional space. We can compute the cosine of the angle between the

two vectors as follows:&lt;/p&gt;

&lt;p&gt;$$ \phi = \text{cos}(\theta)

= \frac{\textbf{u}\cdot\textbf{v}}{||\textbf{u}|| ||\textbf{v}||} $$&lt;/p&gt;

&lt;p&gt;Here &lt;strong&gt;u&lt;/strong&gt; and &lt;strong&gt;v&lt;/strong&gt; are two &lt;em&gt;k-mers&lt;/em&gt; in our complex $k$-dimensional space.

The probability of these two &lt;em&gt;k-mers&lt;/em&gt; being in the same bin, is equal to the

probability that &lt;em&gt;no&lt;/em&gt; random hyperplane cuts through the angle of &lt;strong&gt;u&lt;/strong&gt; and

&lt;strong&gt;v&lt;/strong&gt;, this is visualised in the figure below.&lt;/p&gt;

&lt;p&gt;

&lt;figure &gt;

&lt;img src=&#34;https://lucasvandijk.nl/img/blogimages/kmer-neighbours.png&#34; alt=&#34;

The probability that a hyperplane cuts through the angle of two

similar k-mers is lower than with two very dissimilar k-mers.

&#34; /&gt;

&lt;figcaption&gt;

The probability that a hyperplane cuts through the angle of two

similar k-mers is lower than with two very dissimilar k-mers.

&lt;/figcaption&gt;

&lt;/figure&gt;

&lt;/p&gt;

&lt;p&gt;We end up with the following equation:&lt;/p&gt;

&lt;p&gt;$$ P(h(\textbf{u}) = h(\textbf{v})) = 1 - \frac{\text{acos}(\phi)}{\pi} $$&lt;/p&gt;

&lt;p&gt;Here $h(\textbf{u})$ represent the &lt;em&gt;hash&lt;/em&gt; value of a &lt;em&gt;k-mer&lt;/em&gt; &lt;strong&gt;u&lt;/strong&gt;, or the $n$-bit

integer by generating $n$ random hyperplanes.&lt;/p&gt;

&lt;p&gt;This is the reason why this could be used for &lt;em&gt;approximate&lt;/em&gt; counting of

&lt;em&gt;k-mers&lt;/em&gt;, because some &lt;em&gt;k-mers&lt;/em&gt; may be mistaken for another. By tuning $k$ and

$n$ you can try to minimise the impact on your application. The impact is even

further minimised because similar &lt;em&gt;k-mers&lt;/em&gt; have a higher probability of ending

up in the same bin than dissimilar &lt;em&gt;k-mers&lt;/em&gt;, while reducing the number of bits

you need to store the &lt;em&gt;k-mer&lt;/em&gt;.&lt;/p&gt;

&lt;div class=&#34;footnotes&#34;&gt;

&lt;hr /&gt;

&lt;ol&gt;

&lt;li id=&#34;fn:1&#34;&gt;Zielezinski, A., Vinga, S., Almeida, J. &amp;amp; Karlowski, W. M. &lt;em&gt;Alignment-free sequence comparison: Benefits, applications, and tools.&lt;/em&gt; Genome Biol. 18, 1–17 (2017).

&lt;a class=&#34;footnote-return&#34; href=&#34;#fnref:1&#34;&gt;&lt;sup&gt;^&lt;/sup&gt;&lt;/a&gt;&lt;/li&gt;

&lt;li id=&#34;fn:2&#34;&gt;Cleary, B. et al. &lt;em&gt;Detection of low-abundance bacterial strains in metagenomic datasets by eigengenome partitioning.&lt;/em&gt; Nat. Biotechnol. 33, 1053–1060 (2015).

&lt;a class=&#34;footnote-return&#34; href=&#34;#fnref:2&#34;&gt;&lt;sup&gt;^&lt;/sup&gt;&lt;/a&gt;&lt;/li&gt;

&lt;/ol&gt;

&lt;/div&gt;</description>

</item>

<item>

<title>Visualising the height of the Netherlands</title>

<link>https://lucasvandijk.nl/2018/02/visualising-the-height-of-the-netherlands/</link>

<pubDate>Thu, 01 Feb 2018 00:00:00 +0000</pubDate>

<guid>https://lucasvandijk.nl/2018/02/visualising-the-height-of-the-netherlands/</guid>

<description>&lt;p&gt;On this date 65 years ago, February 1st 1953, the Netherlands experienced its

greatest flood till date, the &lt;a href=&#34;https://en.wikipedia.org/wiki/North_Sea_flood_of_1953&#34; target=&#34;_blank&#34;&gt;North Sea Flood of 1953&lt;/a&gt;. This

is still one of the biggest disasters the Netherlands has ever experienced,

with thousands of casualties and lots of people who lost their homes.&lt;/p&gt;

&lt;p&gt;The Netherlands earns its name from the fact that large parts of the country

lie below sea level. To make sure our country doesn&amp;rsquo;t flood, we have built lots

of barriers, dams and dykes to keep the water out, and we want to prevent

anything like the flood of 1953 from happening ever again.&lt;/p&gt;

&lt;p&gt;In the beginning of this year, several of these barriers and dams were put to

the test when a heavy storm reached the Netherlands which resulted in very high

water levels. Our five biggest dams and barriers needed to be closed at the

same time, a first since their construction. We can ask ourselves the question

whether this will happen more often now that sea levels are rising due to

global warming &lt;sup class=&#34;footnote-ref&#34; id=&#34;fnref:sealevel&#34;&gt;&lt;a href=&#34;#fn:sealevel&#34;&gt;1&lt;/a&gt;&lt;/sup&gt;. A higher base line sea level increases the chance

for even higher water levels when it storms. To get an idea what areas would be

affected the most by a possible flood, we have created a visualisation project

that shows the height of the Netherlands in comparison to the sea level.&lt;/p&gt;

&lt;p&gt;

&lt;figure &gt;

&lt;img src=&#34;https://lucasvandijk.nl/img/blogimages/maeslantkering.jpg&#34; alt=&#34;Maeslantkering&#34; /&gt;

&lt;figcaption&gt;

The Maeslantkering in closed state.

Photo: &lt;a href=&#34;https://beeldbank.rws.nl&#34; target=&#34;_blank&#34;&gt;Rijkswaterstaat&lt;/a&gt;

&lt;/figcaption&gt;

&lt;/figure&gt;

&lt;/p&gt;

&lt;h2 id=&#34;things-to-check-out&#34;&gt;Things to check out&lt;/h2&gt;

&lt;ul&gt;

&lt;li&gt;The Y-axis shows deviation from &amp;ldquo;&lt;a href=&#34;https://en.wikipedia.org/wiki/Amsterdam_Ordnance_Datum&#34; target=&#34;_blank&#34;&gt;Amsterdam Ordnance Datum&lt;/a&gt;&amp;ldquo;&lt;/li&gt;

&lt;li&gt;You can zoom on both the map and the elevation profile (by selecting a region

in the bottom profile).&lt;/li&gt;

&lt;li&gt;If you hover your mouse over the elevation profile, a red dot will show you

the location on the map&lt;/li&gt;

&lt;li&gt;The highest point in our country is the &lt;a href=&#34;https://nl.wikipedia.org/wiki/Vaalserberg&#34; target=&#34;_blank&#34;&gt;Vaalserberg&lt;/a&gt;, with a

whopping 322 m! It&amp;rsquo;s in the most southern part of the Netherlands. Some

people, however, may claim it&amp;rsquo;s not really part of the Netherlands because of

their ugly accents ;)&lt;/li&gt;

&lt;li&gt;If you zoom in on our shorelines, you can see our dykes&lt;/li&gt;

&lt;/ul&gt;

&lt;h2 id=&#34;title&#34;&gt;Elevation profiles of the Netherlands&lt;/h2&gt;

&lt;div id=&#34;chart&#34;&gt;&lt;/div&gt;

&lt;div id=&#34;map&#34;&gt;&lt;/div&gt;

&lt;h2 id=&#34;technical-details&#34;&gt;Technical details&lt;/h2&gt;

&lt;h3 id=&#34;quick-links&#34;&gt;Quick links&lt;/h3&gt;

&lt;ul&gt;

&lt;li&gt;Data source: &lt;a href=&#34;http://ahn.nl&#34; target=&#34;_blank&#34;&gt;http://ahn.nl&lt;/a&gt;

&lt;ul&gt;

&lt;li&gt;We use AHN version 2 and the dataset with a resolution of 5 meter&lt;/li&gt;

&lt;/ul&gt;&lt;/li&gt;

&lt;li&gt;Written in Python and CoffeeScript, and available on Github:

&lt;a href=&#34;https://github.com/lrvdijk/nl-height&#34; target=&#34;_blank&#34;&gt;https://github.com/lrvdijk/nl-height&lt;/a&gt;&lt;/li&gt;

&lt;/ul&gt;

&lt;h3 id=&#34;algorithm-overview&#34;&gt;Algorithm overview&lt;/h3&gt;

&lt;ol&gt;

&lt;li&gt;The publicly available height data is split across multiple &amp;ldquo;raster data&amp;rdquo;

files. The AHN website provides a &amp;ldquo;ahn_units.shp&amp;rdquo; shapefile, that shows

which part of the Netherlands is covered by which &lt;em&gt;raster data&lt;/em&gt; file, and

thus shows which file to download.&lt;/li&gt;

&lt;li&gt;We import this ahn_units shapefile in a PostgreSQL database with the PostGIS

extension enabled. A script downloads and extracts all available raster data

files.&lt;/li&gt;

&lt;li&gt;When all data is downloaded, we import the raster data in our PostgreSQL

database.&lt;/li&gt;

&lt;li&gt;To generate an elevation profile, we draw a line horizontally across the

Netherlands and query which cells in our raster data hit our line. Each cell

contains the height at its midpoint. Our query is based on this blogpost:

&lt;a href=&#34;http://blog.mathieu-leplatre.info/drape-lines-on-a-dem-with-postgis.html&#34; target=&#34;_blank&#34;&gt;http://blog.mathieu-leplatre.info/drape-lines-on-a-dem-with-postgis.html&lt;/a&gt;&lt;/li&gt;

&lt;li&gt;When there&amp;rsquo;s water, no data exists in our database. We apply some

post-processing such that our visualisation does not get weird artefacts.&lt;/li&gt;

&lt;li&gt;We store the data in a JSON file, and use the d3.js library for the

visualisation.&lt;/li&gt;

&lt;/ol&gt;

&lt;div class=&#34;footnotes&#34;&gt;

&lt;hr /&gt;

&lt;ol&gt;

&lt;li id=&#34;fn:sealevel&#34;&gt;&lt;a href=&#34;https://climate.nasa.gov/vital-signs/sea-level/&#34; target=&#34;_blank&#34;&gt;https://climate.nasa.gov/vital-signs/sea-level/&lt;/a&gt;

&lt;a class=&#34;footnote-return&#34; href=&#34;#fnref:sealevel&#34;&gt;&lt;sup&gt;^&lt;/sup&gt;&lt;/a&gt;&lt;/li&gt;

&lt;/ol&gt;

&lt;/div&gt;</description>

</item>

<item>

<title>PHASM: Haplotype-aware de novo genome assembly</title>

<link>https://lucasvandijk.nl/project/phasm/</link>

<pubDate>Tue, 29 Aug 2017 00:00:00 +0000</pubDate>

<guid>https://lucasvandijk.nl/project/phasm/</guid>

<description>

&lt;p&gt;PHASM is a long read &lt;em&gt;de novo&lt;/em&gt; genome assembler that phases variants among

chromosome homologues during the assembly process, and aims to output separate

contigs for each haplotype. The main idea in PHASM is to build bubble chains:

consecutive &amp;ldquo;superbubbles&amp;rdquo; chained together. While most traditional genome

assemblers pop these superbubbles by only keeping the best supported path,

PHASM finds &lt;em&gt;k&lt;/em&gt; paths through this chain of superbubbles that best represent

each haplotype.&lt;/p&gt;

&lt;p&gt;This program has been created as part of my master thesis project. For now, it

has only been tested with error free data.&lt;/p&gt;

&lt;ul&gt;

&lt;li&gt;&lt;strong&gt;Status:&lt;/strong&gt; Finished.&lt;/li&gt;

&lt;li&gt;&lt;strong&gt;License&lt;/strong&gt;: MIT&lt;/li&gt;

&lt;/ul&gt;

&lt;h2 id=&#34;overview&#34;&gt;Overview&lt;/h2&gt;

&lt;h3 id=&#34;pipeline&#34;&gt;Pipeline&lt;/h3&gt;

&lt;p&gt;The PHASM pipeline consists of four main stages:&lt;/p&gt;

&lt;ol&gt;

&lt;li&gt;Overlapper&lt;/li&gt;

&lt;li&gt;Assembly graph construction&lt;/li&gt;

&lt;li&gt;Bubblechain identification&lt;/li&gt;

&lt;li&gt;Phasing&lt;/li&gt;

&lt;/ol&gt;

&lt;p&gt;&lt;img src=&#34;https://lucasvandijk.nl/img/projectimages/phasm-overview.png&#34; alt=&#34;PHASM pipeline overview&#34; /&gt;&lt;/p&gt;

&lt;h3 id=&#34;key-points&#34;&gt;Key Points&lt;/h3&gt;

&lt;ul&gt;

&lt;li&gt;Main idea:

&lt;ul&gt;

&lt;li&gt;Different alleles on homologous chromosomes result in branches in the

assembly graph;&lt;/li&gt;

&lt;li&gt;A bubble or &lt;em&gt;superbubble&lt;/em&gt; is a common motif in a genome assembly graph;&lt;/li&gt;

&lt;li&gt;&lt;strong&gt;Hypothesis&lt;/strong&gt;: each superbubble has a valid path from its source to its

sink that spells a DNA sequence that matches with one of the haplotypes.

Furthermore, we expect that for each haplotype, there exists at least one

path through this superbubble that spells a matching DNA sequence.&lt;/li&gt;

&lt;li&gt;A lot of current genome assemblers pop bubbles in the assembly graph,

thereby removing variation among homologous chromosomes&lt;/li&gt;

&lt;li&gt;Build an assembler that keeps superbubbles and searches for the best set

of paths such that each path represents a haplotype.&lt;/li&gt;

&lt;/ul&gt;&lt;/li&gt;

&lt;li&gt;Main findings:

&lt;ul&gt;

&lt;li&gt;Superbubbles are an incomplete representation of variation among

homologous chromosomes&lt;/li&gt;

&lt;li&gt;Paths through a superbubble are rarely consistent with one of the

haplotypes&lt;/li&gt;

&lt;li&gt;Superbubbles are probably an artefact of approximate overlapping between

reads.&lt;/li&gt;

&lt;/ul&gt;&lt;/li&gt;

&lt;/ul&gt;

&lt;p&gt;Read more in my &lt;a href=&#34;https://lucasvandijk.nl/docs/master-thesis.pdf&#34;&gt;master thesis&lt;/a&gt;.&lt;/p&gt;

&lt;h2 id=&#34;availability&#34;&gt;Availability&lt;/h2&gt;

&lt;p&gt;PHASM is available on &lt;a href=&#34;https://github.com/AbeelLab/phasm&#34; target=&#34;_blank&#34;&gt;Github&lt;/a&gt;.&lt;/p&gt;

</description>

</item>

<item>

<title>Drawing arbitrary shapes with OpenGL points</title>

<link>https://lucasvandijk.nl/2015/06/drawing-arbitrary-shapes-with-opengl-points/</link>

<pubDate>Sat, 13 Jun 2015 00:00:00 +0000</pubDate>

<guid>https://lucasvandijk.nl/2015/06/drawing-arbitrary-shapes-with-opengl-points/</guid>

<description>&lt;p&gt;Part of my Google Summer of Code project involves porting several arrow

heads from &lt;a href=&#34;https://github.com/glumpy/glumpy&#34; target=&#34;_blank&#34;&gt;Glumpy&lt;/a&gt; to

&lt;a href=&#34;https://vispy.org&#34; target=&#34;_blank&#34;&gt;Vispy&lt;/a&gt;. I also want to make a slight change to them:

the arrow heads in Glumpy include an arrow body, I want to remove that

to make sure you can put an arrow head on every type of line you want.&lt;/p&gt;

&lt;p&gt;Making a change like that requires that you understand how those shapes

are drawn. And for someone without a background in computer graphics this took some thorough

investigation of the code and the techniques used. This article is aimed

at people like me: good enough programming skills and linear algebra

knowledge, but almost no former experience with OpenGL or computer

graphics in general.&lt;/p&gt;

&lt;h2 id=&#34;implicit-surfaces&#34;&gt;Implicit surfaces&lt;/h2&gt;

&lt;p&gt;The basic principle behind drawing these 2D shapes is called implicit

surfaces, and it relies on a arbitrary shape function which for a given

point in your image returns the distance to your shape surface or

boundary from that given point. This is visualized in fig. 1.&lt;/p&gt;

&lt;p&gt;

&lt;figure &gt;

&lt;img src=&#34;https://lucasvandijk.nl/img/blogimages/opengl-points/implicit-surfaces.png&#34; alt=&#34;

Distances to the boundary of a circle

&#34; /&gt;

&lt;figcaption&gt;

Distances to the boundary of a circle

&lt;/figcaption&gt;

&lt;/figure&gt;

&lt;/p&gt;

&lt;p&gt;The distances any shape function should return are highlighted in red.

To actually be able to draw a shape, we need to distinguish whether a

point lies inside the shape or not. We make the arbitrary decision that

a negative distance lies inside a shape, and a positive distance means

that the point lies outside the shape.&lt;/p&gt;

&lt;h2 id=&#34;distance-functions-for-a-few-basic-shapes&#34;&gt;Distance functions for a few basic shapes&lt;/h2&gt;

&lt;p&gt;The distance functions defined below have one requirement: the center

point of the shape has the coordinate (0, 0).&lt;/p&gt;

&lt;h3 id=&#34;circle&#34;&gt;Circle&lt;/h3&gt;

&lt;p&gt;For a circle these distances are easily calculated:&lt;/p&gt;

&lt;p&gt;$$ d(\textbf{x}) = ||\textbf{x}|| − r $$&lt;/p&gt;

&lt;p&gt;Where:&lt;/p&gt;

&lt;ul&gt;

&lt;li&gt;$\textbf{x}$: Vector to the point in consideration.&lt;/li&gt;

&lt;li&gt;$r$: The radius of the circle.&lt;/li&gt;

&lt;/ul&gt;

&lt;p&gt;If the point lies within the circle, the length of the vector towards

that point is smaller than the radius, making the distance automatically

negative.&lt;/p&gt;

&lt;h3 id=&#34;square&#34;&gt;Square&lt;/h3&gt;

&lt;p&gt;Consider fig. 2.&lt;/p&gt;

&lt;p&gt;

&lt;figure &gt;

&lt;img src=&#34;https://lucasvandijk.nl/img/blogimages/opengl-points/distance-square.png&#34; alt=&#34;

Visualizing the distances to the border of a square

&#34; /&gt;

&lt;figcaption&gt;

Visualizing the distances to the border of a square

&lt;/figcaption&gt;

&lt;/figure&gt;

&lt;/p&gt;

&lt;p&gt;A square is a nice symmetric figure, so we can take the absolute values

of the coordinates, and then we only need to consider the smaller

highlighted square (light blue).&lt;/p&gt;

&lt;p&gt;The distance to the boundary of the square is then:&lt;/p&gt;

&lt;p&gt;$$ d(\textbf{x}) = \text{max}(|x_1|, |x_2|) - \frac{s}{2}$$&lt;/p&gt;

&lt;p&gt;Where:&lt;/p&gt;

&lt;ul&gt;

&lt;li&gt;$\textbf{x}$: Vector towards the point in consideration&lt;/li&gt;

&lt;li&gt;$|x_1|, |x_2|$ are the absolute values of the first and second element of the

vector (the &lt;em&gt;x&lt;/em&gt; and &lt;em&gt;y&lt;/em&gt; coordinates).&lt;/li&gt;

&lt;li&gt;$s$: The size of the square.&lt;/li&gt;

&lt;/ul&gt;

&lt;p&gt;Using the max function we sort of select to which boundary the distance

will be calculated. Then we substract the size of the smaller square

(highlighted with light blue). The resulting distance is then negative

if the point lies within the square, and positive otherwise.&lt;/p&gt;

&lt;h2 id=&#34;combining-shapes&#34;&gt;Combining shapes&lt;/h2&gt;

&lt;p&gt;To make more complex shapes, it is often useful to combine multiple

simple shapes. To do this, we introduce some basic operations on the

returned distances of a simple shape &lt;sup class=&#34;footnote-ref&#34; id=&#34;fnref:1&#34;&gt;&lt;a href=&#34;#fn:1&#34;&gt;1&lt;/a&gt;&lt;/sup&gt;.&lt;/p&gt;

&lt;h3 id=&#34;inverse&#34;&gt;Inverse&lt;/h3&gt;

&lt;p&gt;We made the arbitrary decision to say that the distance is negative when

a point lies within the shape. To get the inverse of a shape we simply

need to negate each distance value.&lt;/p&gt;

&lt;p&gt;$$\forall x, y: \neg S(x, y) = -S(x, y)$$&lt;/p&gt;

&lt;h3 id=&#34;union&#34;&gt;Union&lt;/h3&gt;

&lt;p&gt;The union of two shapes can be retrieved by using the min function on

both distance functions.&lt;/p&gt;

&lt;p&gt;$$ \forall x, y : U(x, y) = \min(S_1(x, y), S_2(x, y)) $$&lt;/p&gt;

&lt;p&gt;Remember that the distance value is negative when the point lies in the

shape. The lowest value will win here, so if one of those distances is

negative (the point belongs to at least one shape), it will return the

negative value. Thus combining both shapes to a single one.&lt;/p&gt;

&lt;h3 id=&#34;difference&#34;&gt;Difference&lt;/h3&gt;

&lt;p&gt;The difference of two shapes contains all points in $S_1$ excluding the

points in $S_2$. This is defined as follows:&lt;/p&gt;

&lt;p&gt;$$\forall x, y : D(x, y) = \max(S_1(x, y), -S_2(x, y))$$&lt;/p&gt;

&lt;p&gt;Consider the example where $S_1(x, y)=−2$ and

$S_2(x, y)=−1$. In short, the current point $(x, y)$ belongs to

both $S_1$ and $S_2$. Using the above function for the difference, the value

from $S_2$ gets negated: $−S_2(x, y) = 1$. Due to the max function,

this value will win (it&amp;rsquo;s higher than -2), and therefore, it will &lt;em&gt;not&lt;/em&gt;

be part of the new shape. Which is precisely what we want, because we

want all points that are part of $S_1$ but not part of $S_2$.&lt;/p&gt;

&lt;h3 id=&#34;intersection&#34;&gt;Intersection&lt;/h3&gt;

&lt;p&gt;The intersection of two shapes contains all points that are both part of

$S_1$ and $S_2$. It is defined as follows:&lt;/p&gt;

&lt;p&gt;$$\forall x, y: I(x, y) = \max(S_1(x, y), S_2(x, y))$$&lt;/p&gt;

&lt;p&gt;A point will be part of the new shape if and only if both distances are

negative. If one distance is positive, the max function will return this

value, and a positive value means not part of the shape. This results in

a shape which includes only points that are both part of $S_1$ and $S_2$.&lt;/p&gt;

&lt;h2 id=&#34;opengl-implementation&#34;&gt;OpenGL implementation&lt;/h2&gt;

&lt;p&gt;So how do we translate these principles to some working code? Let us

first introduce some basic OpenGL concepts before throwing the shader

code at you.&lt;/p&gt;

&lt;h3 id=&#34;shaders-and-drawing-modes&#34;&gt;Shaders and drawing modes&lt;/h3&gt;

&lt;p&gt;I will not go too deep in the basics of OpenGL, but a modern OpenGL

pipeline consists of multiple &lt;em&gt;shaders&lt;/em&gt;, small programs you can write

yourself. At the very minimum you&amp;rsquo;ll need a &lt;em&gt;vertex shader&lt;/em&gt; and a

&lt;em&gt;fragment shader&lt;/em&gt;, which determine where the main &amp;ldquo;drawing points&amp;rdquo; will

be positioned and the color of the individual pixels respectively. The

pipeline is illustrated in fig. 3.&lt;/p&gt;

&lt;p&gt;

&lt;figure &gt;

&lt;img src=&#34;https://lucasvandijk.nl/img/blogimages/opengl-points/gl-pipeline.png&#34; alt=&#34;

OpenGL pipeline illustrated.

&#34; /&gt;

&lt;figcaption&gt;

OpenGL pipeline illustrated.

&lt;/figcaption&gt;

&lt;/figure&gt;

&lt;/p&gt;

&lt;p&gt;You can define your own attributes for each vertex, for example the

position, colour, or orientation.&lt;/p&gt;

&lt;p&gt;OpenGL has several modes for generating the actual primitives. These are

illustrated in fig. 4.&lt;/p&gt;

&lt;p&gt;

&lt;figure &gt;

&lt;img src=&#34;https://lucasvandijk.nl/img/blogimages/opengl-points/gl-primitives.png&#34; alt=&#34;

OpenGL primitive generation modes.

&#34; /&gt;

&lt;figcaption&gt;

OpenGL primitive generation modes.

&lt;/figcaption&gt;

&lt;/figure&gt;

&lt;/p&gt;

&lt;p&gt;For a more in depth OpenGL introduction, I would recommend &lt;a href=&#34;http://www.labri.fr/perso/nrougier/teaching/opengl/&#34; target=&#34;_blank&#34;&gt;this

tutorial&lt;/a&gt;, &lt;a href=&#34;http://antongerdelan.net/opengl/&#34; target=&#34;_blank&#34;&gt;Anton&amp;rsquo;s

OpenGL tutorials&lt;/a&gt;, or

&lt;a href=&#34;http://opengl-tutorial.org&#34; target=&#34;_blank&#34;&gt;opengl-tutorial.org&lt;/a&gt;.&lt;/p&gt;

&lt;p&gt;&lt;small&gt;Above pictures courtesy of N. Rougier&lt;sup class=&#34;footnote-ref&#34; id=&#34;fnref:2&#34;&gt;&lt;a href=&#34;#fn:2&#34;&gt;2&lt;/a&gt;&lt;/sup&gt;&lt;/small&gt;&lt;/p&gt;

&lt;h3 id=&#34;general-2d-shape-shaders&#34;&gt;General 2D shape shaders&lt;/h3&gt;

&lt;p&gt;For the 2D shapes we want to draw, we&amp;rsquo;re going to use the points drawing

mode. OpenGL allows you to specify the point size, and the fragment

shader will be called for each pixel in the point.&lt;/p&gt;

&lt;p&gt;Let&amp;rsquo;s check the vertex shader were we position our vertices, and

configure the point size.&lt;/p&gt;

&lt;h4 id=&#34;vertex-shader&#34;&gt;Vertex shader&lt;/h4&gt;

&lt;p&gt;Listing 1: The vertex shader code for our 2D shapes&lt;/p&gt;

&lt;pre&gt;&lt;code class=&#34;language-glsl&#34;&gt;// Uniforms

// ------------------------------------

uniform float antialias;

uniform mat4 ortho;

// Attributes

// ------------------------------------

attribute vec2 position;

attribute float size;

attribute vec4 color;

attribute float rotation;

attribute float linewidth;

// Varyings

// ------------------------------------

varying float v_size;

varying vec4 v_color;

varying vec2 v_rotation;

varying float v_antialias;

varying float v_linewidth;

// Main

// ------------------------------------

void main (void)

{

v_size = size;

v_linewidth = linewidth;

v_antialias = antialias;

v_color = color;

v_rotation = vec2(cos(rotation), sin(rotation));

gl_Position = ortho * vec4(position, 0, 1);

gl_PointSize = M_SQRT2 * size + 2.0 * (linewidth + 1.5*antialias);

}

&lt;/code&gt;&lt;/pre&gt;

&lt;p&gt;We first define some &lt;em&gt;uniforms&lt;/em&gt;, &lt;em&gt;attributes&lt;/em&gt;, and &lt;em&gt;varyings&lt;/em&gt;. Uniforms

are variables which are the same for each vertex. Attributes are

variables defined per vertex, and with varyings we can pass some data to

the next steps in the pipeline (for example, the fragment shader).&lt;/p&gt;

&lt;p&gt;Each vertex has a position, and this is where our 2D shape will be

drawn. The matrix &lt;code&gt;ortho&lt;/code&gt; is used for the proper projection of the

vertex to your screen. We will not explain this in-depth, but if you

want to know more please refer to &lt;a href=&#34;http://www.opengl-tutorial.org/beginners-tutorials/tutorial-3-matrices/&#34; target=&#34;_blank&#34;&gt;this

tutorial&lt;/a&gt;

on opengl-tutorial.org.&lt;/p&gt;

&lt;p&gt;Our shapes are larger than one pixel, so we need to change

&lt;code&gt;gl_PointSize&lt;/code&gt;. Our shapes also have a size attached to them, but for

the point size we scale this size with $\sqrt{2}$ (ignore the extra

size for antialias en linewidth for now). We do this because our shapes

can be rotated. To fit a rotated square of size &lt;em&gt;x&lt;/em&gt; in another square,

we need a bigger square of size $x\sqrt{2}$ (I hope you guys remember

Pythagoras). This is visualized in fig. 5.&lt;/p&gt;

&lt;p&gt;

&lt;figure &gt;

&lt;img src=&#34;https://lucasvandijk.nl/img/blogimages/opengl-points/point-rotation.png&#34; alt=&#34;

Rotation of a square

&#34; /&gt;

&lt;figcaption&gt;

Rotation of a square

&lt;/figcaption&gt;

&lt;/figure&gt;

&lt;/p&gt;

&lt;p&gt;Furthermore, we pass along some variables to the next steps in the

pipeline (size, line width, antialias, color, rotation). Note we create

a direction vector for the rotation from the given rotation in radians.&lt;/p&gt;

&lt;h4 id=&#34;fragment-shader&#34;&gt;Fragment shader&lt;/h4&gt;

&lt;p&gt;Listing 2: The fragment shader code&lt;/p&gt;

&lt;pre&gt;&lt;code class=&#34;language-glsl&#34;&gt;// Varyings

// ------------------------------------

varying float v_size;

varying vec4 v_color;

varying vec2 v_rotation;

varying float v_antialias;

varying float v_linewidth;

// Main

// ------------------------------------

void main()

{

vec2 P = gl_PointCoord.xy - vec2(0.5,0.5);

P = vec2(v_rotation.x*P.x - v_rotation.y*P.y,

v_rotation.y*P.x + v_rotation.x*P.y) * v_size;

float size = v_size/M_SQRT2;

float distance = shape_func(P, size);

gl_FragColor = filled(distance, v_linewidth, v_antialias, v_color);

}

&lt;/code&gt;&lt;/pre&gt;

&lt;p&gt;Note we have the same varying definitions as in the vertex shader. These

contain values as passed from the vertex shader. Also note the usage of

the built in variable &lt;code&gt;gl_PointCoord&lt;/code&gt;. We specified in the vertex shader

the size of our point in pixels, and for each pixel in this point the

fragment shader gets called. The &lt;code&gt;gl_PointCoord&lt;/code&gt; contains the

coordinates &lt;strong&gt;inside the point&lt;/strong&gt;. The &lt;em&gt;x&lt;/em&gt; and &lt;em&gt;y&lt;/em&gt; attributes from

&lt;code&gt;gl_PointCoord&lt;/code&gt; range from 0.0 to 1.0, where (0, 0) is the bottom left

corner of the point, and (1, 1) is the top right corner of the point.&lt;/p&gt;

&lt;p&gt;There are several operations applied to these coordinates:&lt;/p&gt;

&lt;p&gt;&lt;strong&gt;Step 1.&lt;/strong&gt; First we substract $\begin{bmatrix}0.5 \\ 0.5 \end{bmatrix}$.

This makes sure the origin is in the center of the point, because

the distance functions we defined earlier in this article

require that.&lt;/p&gt;

&lt;p&gt;&lt;strong&gt;Step 2.&lt;/strong&gt; Next we apply a rotation transformation to the point.&lt;/p&gt;

&lt;p&gt;Remember the transformation matrix is:&lt;/p&gt;

&lt;p&gt;$$ \begin{bmatrix} \text{cos}(\theta) &amp;amp; \text{-sin}(\theta) \\

\text{sin}(\theta) &amp;amp; \text{cos}(\theta) \end{bmatrix}$$&lt;/p&gt;

&lt;p&gt;If you look closely at the code you see that the vector &lt;code&gt;P&lt;/code&gt; gets

multiplied with this matrix.&lt;/p&gt;

&lt;p&gt;$$\begin{bmatrix}nx \\ ny\end{bmatrix} = \begin{bmatrix}x \\ y

\end{bmatrix} \begin{bmatrix} \mathrm{cos}(\theta) &amp;amp; \mathrm{-sin}(\theta) \\

\mathrm{sin}(\theta) &amp;amp; \mathrm{cos}(\theta) \end{bmatrix} = \begin{bmatrix} x

\mathrm{cos}(\theta) - y \mathrm{sin}(\theta) \\ x

\mathrm{sin}(\theta) + y \mathrm{cos}(\theta) \end{bmatrix} $$&lt;/p&gt;

&lt;p&gt;&lt;strong&gt;Step 3.&lt;/strong&gt; We also scale the coordinates with &lt;code&gt;v_size&lt;/code&gt;.&lt;/p&gt;

&lt;p&gt;These transformations are visualized in fig. 6.&lt;/p&gt;

&lt;p&gt;

&lt;figure &gt;

&lt;img src=&#34;https://lucasvandijk.nl/img/blogimages/opengl-points/transformations.png&#34; alt=&#34;

Transformation of the coordinates within a point

&#34; /&gt;

&lt;figcaption&gt;

Transformation of the coordinates within a point

&lt;/figcaption&gt;

&lt;/figure&gt;

&lt;/p&gt;

&lt;p&gt;The green region in the bottom of fig. 6 is our canvas for drawing our

shape. This is done by &lt;code&gt;shape_func()&lt;/code&gt;, any function that returns the

distance to the boundary of a shape as explained earlier in this

article.&lt;/p&gt;

&lt;p&gt;The &lt;code&gt;filled()&lt;/code&gt; function determines the color for the current pixel

determined by the returned distance of &lt;code&gt;shape_func()&lt;/code&gt;. Simply put: if

the returned distance is negative, it returns a color (because it&amp;rsquo;s part

of the shape), otherwise it makes the current pixel transparent. It also

applies some anti-aliasing techniques which I don&amp;rsquo;t know the details of,

so we will not cover this in-depth.&lt;/p&gt;

&lt;h2 id=&#34;example-curved-arrows&#34;&gt;Example: curved arrows&lt;/h2&gt;

&lt;p&gt;To conclude this article we will look into the distance function of a

curved arrow head. A curved arrow head can be constructed using the

inverse of the union of three circles. This is visualized in fig. 7, and

the accompanying shader code can be seen in lst. 3.&lt;/p&gt;

&lt;p&gt;

&lt;figure &gt;

&lt;img src=&#34;https://lucasvandijk.nl/img/blogimages/opengl-points/curved-arrows.png&#34; alt=&#34;

Construction of a curved arrow head

&#34; /&gt;

&lt;figcaption&gt;

Construction of a curved arrow head

&lt;/figcaption&gt;

&lt;/figure&gt;

&lt;/p&gt;

&lt;p&gt;Listing 3: GLSL function to the distance of an arrow&lt;/p&gt;

&lt;pre&gt;&lt;code class=&#34;language-glsl&#34;&gt;/**

* Computes the signed distance to an curved arrow

*

* Parameters:

* -----------

* texcoord : Point to compute distance to

* size : size of the arrow head in pixels

*

* Return:

* -------

* Signed distance to the arrow

*

*/

float arrow_curved(vec2 texcoord, float size)

{

vec2 start = -vec2(size/2, 0.0);

vec2 end = +vec2(size/2, 0.0);

float height = 0.5;

vec2 p1 = start + size*vec2(0, -height);

vec2 p2 = start + size*vec2(0, +height);

vec2 p3 = end;

// Head : 3 circles

vec2 c1 = circle_from_2_points(p1, p3, 6.0*size).zw;

float d1 = length(texcoord - c1) - 6*size;

vec2 c2 = circle_from_2_points(p2, p3, 6.0*size).xy;

float d2 = length(texcoord - c2) - 6*size;

vec2 c3 = circle_from_2_points(p1, p2, 3.0*size).xy;

float d3 = length(texcoord - c3) - 3*size;

return -min(d3, min(d1,d2));

}

&lt;/code&gt;&lt;/pre&gt;

&lt;p&gt;We first define the arrow corner points as &lt;code&gt;p1&lt;/code&gt;, &lt;code&gt;p2&lt;/code&gt; and &lt;code&gt;p3&lt;/code&gt;. We use a

helper function which calculates the center point of a circle through

two points with a given radius &lt;code&gt;r&lt;/code&gt;.&lt;/p&gt;

&lt;p&gt;The distance to these circles are then easily calculated, and we use the

min function to get the union of these three circles. Our arrow head is

exactly the area &lt;em&gt;not&lt;/em&gt; covered by these circles (see fig. 7), and

therefore we return the inverse of this union.&lt;/p&gt;

&lt;h3 id=&#34;references&#34;&gt;References&lt;/h3&gt;

&lt;div class=&#34;footnotes&#34;&gt;

&lt;hr /&gt;

&lt;ol&gt;

&lt;li id=&#34;fn:1&#34;&gt;&lt;p&gt;N. P. Rougier, “Antialiased 2D grid, marker, and arrow shaders,”

&lt;em&gt;Journal of Computer Graphics Techniques (JCGT)&lt;/em&gt;, vol. 3, no. 4, pp.

1–52, November 2014 [Online]. Available:

&lt;a href=&#34;http://jcgt.org/published/0003/04/01/&#34; target=&#34;_blank&#34;&gt;http://jcgt.org/published/0003/04/01/&lt;/a&gt;&lt;/p&gt;

&lt;a class=&#34;footnote-return&#34; href=&#34;#fnref:1&#34;&gt;&lt;sup&gt;^&lt;/sup&gt;&lt;/a&gt;&lt;/li&gt;

&lt;li id=&#34;fn:2&#34;&gt;&lt;p&gt;“Modern openGL tutorial with python.” [Online]. Available: