\documentclass[swedish, a4paper]{article}
\usepackage[latin1]{inputenc}
%\usepackage{isolatin1}
\usepackage[T1]{fontenc}
\usepackage{babel}
\usepackage{amstex}
\title{Raytracing med simplexalgoritmen}
\author{Niels Möller <nisse@@lysator.liu.se>}

\frenchspacing
%\raggedright
\setcounter{secnumdepth}{0}

\newcommand{\VEC}[1]{\overline{#1}}

\begin{document}
\maketitle

\vspace{\stretch{1}}
\tableofcontents
\vspace{\stretch{2}}
\newpage

\section{Bakgrund}
Fredrik Hübinette <profezzorn@@lysator.liu.se> och jag pratade för ett
tag sen om raytracing. Hubbe hade en idé om att göra en förprocessor
till en raytracer/doomklon; tanken är att utgående från en lista av
trianglar dela upp rummet i ett antal konvexa polyedrar, och tabellera
dessa. För att göra den uppdelningen behöver man bland annat lösa
problemet att avgöra om en triangel skär en polyeder, avgränsad av
plan. Och just den sortens konvexa polyedrar är simplexalgoritmen bra på
att hantera.

Plan representeras av sin normalvektor $\VEC{n}$ och sitt avstånd $c$
till origo:
\begin{equation}
P = \{ \VEC{r} = (x, y, z) : \VEC{r} \cdot \VEC{n} = c \}
\end{equation}

Ett halvrum representeras på samma sätt, men med likheten ersatt av en
olikhet. En konvex polyeder ges av snittet mellan att antal sådana
halvrum. För att det ska vara enkelt att använda simplexalgoritmen
krävs dessutom att $x,y,z,c \geq 0$. Att $|n| = 1$ är däremot inte
viktigt.

En triangel representeras av det plan som den ligger i, tillsammans
med ytterligare tre plan som avgränsar triangeln. Med andra ord, som
snittet mellan ett plan och tre halvrum.

Givet ett antal trianglar och ett antal konvexa polyedrar, så kan man
göra så här (i pseudopythonkod)
\begin{verbatim}
for t in triangles:
   for p in polyeders:
      if intersectp(p, t):
         p1, p2 = split(p, t)
         polyeders.replace(p, p1, p2)
\end{verbatim}
Här betyder \texttt{split} att låta det plan som triangeln ligger i
skära polyderern i två delar, vilket är rätt enkelt. Det kluriga är
\texttt{intersectp}, att avgöra om en triangel skär en polyeder eller
inte. Den här artikeln försöker reda ut hur detta problem kan lösas
med hjälp av simplexalgoritmen.


\section{Simplexanpassad problemformulering}

Vi vill undersöka snittet mellan en triangel och en konvex
polyeder. Snittet ges av lösningsmängden till
\begin{eqnarray}
  \label{eq:triangelns plan} \VEC{r} \cdot \VEC{n} & = & c \\
  \label{eq:halvrum<} \VEC{r} \cdot \VEC{a_i} & \leq & \alpha_i
  \text{ , för $1 \leq i \leq k$} \\
  \label{eq:halvrum>} \VEC{r} \cdot \VEC{b_i} & \geq & \beta_i
  \text{ , för $1 \leq i \leq l$}
\end{eqnarray}
Här ska samtliga $c, \alpha_i, \beta_i \geq 0$. Ekvation
(\ref{eq:triangelns plan}) beskriver planet som triangeln ligger i,
medan olikheterna beskriver både polyedern och avgränsningen av
triangeln.

Första steget är att införa så kallade \emph{slackvariabler}, för att
kunna ersätta olikheterna med likheter. Då ersätts ekvationerna
(\ref{eq:halvrum<}) och (\ref{eq:halvrum>}) av
\begin{eqnarray}
  \label{eq:slack+} \VEC{r} \cdot \VEC{a_i} + s_i & = & \alpha_i
  \text{ , för $1 \leq i \leq k$} \\
  \label{eq:slack-} \VEC{r} \cdot \VEC{b_i} - r_i & = & \beta_i
  \text{ , för $1 \leq i \leq l$}
\end{eqnarray}
Även $r_i, s_i \geq 0$.

Vi har nu $m = 1 + k + l$ ekvationer, i de $3 + k + l$ variablerna
$\{x, y, z, s_i, r_i\}$. I simplexalgoritmen hanterar man
\emph{tillåtna baslösningar}, vilket betyder att man sätter alla
varibler utom $m$ stycken till 0 (så kallade \emph{icke-basvariabler})
och löser ut de övriga $m$ \emph{basvariablerna} ur ekvationssystemet.
För att resultatet ska ligga i det tillåtna området krävs dels att det
ekvationssystem man får när man stryker icke-basvariablerna är
inverterbart, dels att elementen i lösningsvektorn av basvariabler är
icke-negativa.

I vårt ekvationssystem är det mest praktiskt att starta med $x=y=z=0$
och $s_i = \alpha_i$. Det vill säga, $s_i$:na är basvariabler, och $x,
y, x$ är icke-basvariabler. För att få tillräckligt många behändiga
basvariabler inför vi så kallade \emph{artificiella variabler} i
ekvationerna (\ref{eq:triangelns plan}) och (\ref{eq:slack-}). Då får
vi ett nytt ekvationssystem,
\begin{eqnarray}
\VEC{r} \cdot \VEC{n} + w_0  & = & c \\
\VEC{r} \cdot \VEC{a_i} + s_i & = & \alpha_i , 1 \leq i \leq k \\
\VEC{r} \cdot \VEC{b_i} - r_i + w_i & = & \beta_i, 1 \leq i \leq l \\
x, y, z, s_i, r_i, w_i & \geq & 0
\end{eqnarray}
En lösning till detta ekvationssystem ger en punkt i det snitt vi är
intresserade av, om och endast om alla $w_i = 0$. Därför söker vi
minimum av 
\begin{equation}
\label{eq:ml} f = \sum_0^l w_i
\end{equation}
Triangeln och polyedern skär varandra om och endast om detta minimum
är 0.

Vi har nu $m = 1 + k + l$ ekvationer i $1 + k + 2l$ variabler, och för
detta nya problem kan vi enkelt välja en startpunkt: välj som
basvariabler $s_i = \alpha_i$ för $1 \leq i \leq k$, $w_0 = c$, $w_i =
\beta_i$ för $1 \leq i \leq l$

Innan vi startar simplexalgoritmen måste vi eliminera basvariablerna
ur målfunktionen. Om vi löser ut $w_i$:na, sätter in i (\ref{eq:ml}),
och flyttar över konstanterna, får vi 
\begin{equation}
\label{eq:kanoniskt ml}
f + \VEC{r} \cdot (\VEC{n} + \sum_1^l \VEC{b_i}) - \sum_1^l r_i = c +
\sum_1^l \beta_i
\end{equation}

Hela tiden i simplexalgoritmen har man delat upp variablerna i
$m$ basvariabler och $n-m$ icke-basvariabler, vilket svarar mot ett
visst hörn i det tillåtna området. Och man eliminerar basvariablerna
ur målfunktionen.

För att minska värdet på $f$ i ekvation (\ref{eq:kanoniskt ml}), så tar
vi någon av de variabler vars koefficient är positiv, och ökar
den. Detta är vår \emph{inkommande} basvariabel. Vi
kan öka den ända tills någon av basvariablerna (som inte ingår i
målfunktionen, utan bara i bivillkoren) blir 0. De. Vi tar alltså och
byter ut en av basvariablerna, på ett sådant sätt att målfunktionens
värde minskar. På samma gång så flyttar vi oss längs en kant i det
tillåtna området från ett hörn till ett grannhörn. För detaljer, se
nästa avsnitt.

Simplexalgoritmen kan stoppa på två sätt
\begin{itemize}
\item Ingen inkommande basvariabel: Om ingen av icke-basvariablerna
  har en positiv koefficient i uttrycket för målfunktionen, så betyder
  att det inte går att minska målfunktionen ytterligare. I vårt fall
  så skär triangeln och polyedern varandra om och endast om
  målfunktionens värde är 0. Men eftersom man i praktiken använder
  flyttal, så får man i stället använda villkoret $f < \epsilon$.
\item Ingen utgående basvariabel: Om den nya basvariabeln kan ökas
  obegränsat utan att nån av de andra variablerna blir 0, så är
  målfunktionen nedåt obegränsad. I vårt fall kan det egentligen inte
  inträffa; målfunktionen $\sum w_i$ kan aldrig bli negativ. Men
  konstiga avrundningsfel skulle möjligen kunna stoppa simplexalgoritmen
  här i alla fall.
\end{itemize}


\section{Simplexalgoritmen}

För att lösa ett minimeringsproblem med $n$ variabler och $m$
bivillkor, så använder vi följande data: en $n$-vektor $c$,
en $m \times n$-matris $A$, ett högerled $b$ som är en vektor av $m$
icke-negativa tal, och en indexvektor $B$ av $m$ index som beskriver
de aktuella basvariablerna. Dessutom behöver vi ett ''högerled'' till
målfunktionen, som vi kallar $b_0$ (detta tal kan i allmänhet vara
negativt, men inte med den målfunktion vi använder). Problemet kan då
formuleras
\begin{eqnarray}
  \label{eq:mlrad}\min f \text{ då } f + c^T x & = & b_0 \\
  Ax & = & b \\ 
  x_i & \geq & 0 \text{ , för  $1 \leq i \leq n$}
\end{eqnarray}

När vi arbetar med problemet vill vi hela tiden ha det på
\emph{kanonisk form}, vilket betyder att basvariablerna ska vara
eliminerade ur (\ref{eq:mlrad}), och att motsvarande kolonner i $A$
ska bilda en enhetsmatris. Med andra ord ska $A_{kB_j} =
\delta_{kj}$. Man behöver förstås inte spara dessa koloner i A, och i
vilket fall så måste man räkna som om elementen där är \emph{exakt} 0
och 1.

Simplexalgoritmen går ut på att genomföra följande tre steg ända tills
nåt av stoppvillkoren är uppfyllt.

\begin{itemize}
\item Välj ny inkommande basvariabel.
\item Välj utgående basvariabel.
\item Gör radoperationer i $A$, $b$ och $c$ tills problemet är på kanonisk
form (för den nya uppsättningen basvariabler).
\end{itemize}

Ny inkommande basvariabler måste väljas bland de som har positiva
koefficienter i målfunktionen, det vill säga i $c$. Men vilken ska man
ta om det finns flera? En variant är att helt enkelt ta den som är
störst. Men då kan man, om man har riktig otur, råka utför att
algoritmen börjar cykla. Vill man försäkra sig mot cykling så kan man
i stället ta den av variablerna som har lägst index (kallas
\emph{Blands regel}).

Om vi har valt $x_r$ som inkommande basvariabel, så väljer vi som
utgående basvariabel den som först blir 0 när vi ökar $x_r$. Det är
bara de basvariabler som har $A_{kr} > 0$ som kommer att minska. Och
av dessa så är det den basvariabel $x_{B_k}$ för vilken $b_k/A_{kr}$
är minst som först blir 0. Med andra ord, bestäm det $k$ som uppfyller
$A_{kr} > 0$ och minimerar $b_{k}/A_{kr}$. Då blir $x_{B_k}$ utgående
basvariabel, och vi kommer ihåg det genom att sätta $B_k = r$.

Nu har vi ändrat uppsättningen av basvariabler, genom att sätta $B_k =
r$. För att få problemet på kanonisk form så eliminerar vi element ur
$A$ på ungefär samma sätt som vid vanlig gausselimination. Först delar
vi rad $k$ i $A$-matrisen med $A_{kr}$, så att $A_{kr} = 1$. Sen
subtraherar vi multipler av denna rad från övriga rader i $A$, så att
resten av elementen i kolonn $r$ blir 0. På samma sätt elimineras
$c_r$ ur målfunktionsraden. Och så måste vi förstås göra likadant med
högerleden.

Om $A, B, c, b, b_0$ är de ursprungliga värdena, och $A', B', c', b',
b'_0$ betecknar de nya värdena, så kan en iteration beskrivs så här:

\begin{enumerate}
\item Om alla $c_i \leq 0$, stoppa.
\item Välj ett $r$, med $c_r > 0$ (med eller utan Blands regel).
\item Om alla $A_{jr} \leq 0$, stoppa. 
\item Tag minimum av $\{ b_j / A_{jr} | A_{jr} > 0\}$, och låt $k$ vara det
lägsta index som ger minimum.
\item Sätt $B'_j = B_j$ för $j \neq k$, $B'_k = r$.
\item Sätt $A'_k = A_k / A_{kr}$ och $b'_k =
b_k/A_{kr}$. ($A_k$ betecknar här rad $k$ i $A$).
\item För $j \neq k$, sätt $A'_j = A_j - A_{jr} A_k$ och $b'_j = b_j -
A_{jr} b_k$.
\item Sätt $c' = c - c_r A_k$ och $b'_0 = b_0 - c_r b_k$
\end{enumerate}

Det antal steg som behövs för att hitta minimum är oftast mindre än
$3m$, där $m$ är antalet bivillkor.


\section{Exempel}

Vi vill undersöka om polyedern
\begin{eqnarray} 
x + y + z & \leq & 2 \\
x, y, z & \geq & 0
\end{eqnarray}
skär triangeln 
\begin{eqnarray} 
  x - z & = & 1 \\
  y & \geq & 1 \\
  x & \leq & 3 \\
  -x + y & \leq & 1
\end{eqnarray}
Inför vi slackvariabler i olikheterna får vi ekvationssystemet
\begin{eqnarray}
x + y + z + s_1 & = & 2 \\
\label{eq:ex plan}x - z & = & 1 \\
\label{eq:y>}y - r & = & 1 \\
x + s_2 & = & 3 \\
-x + y + s_3 & = & 1 \\
x, y, z, s_1, s_2, s_3, r & \geq & 0
\end{eqnarray}

Vi saknar basvariabler i ekvationerna (\ref{eq:ex plan}) och
(\ref{eq:y>}), och behöver därför införa två artificiella variabler
$w_1, w_2$. 
\begin{equation}
  \begin{array}{c|rrrrrrrrr|c}
        &  x &  y &  z & s_1 & s_2 & s_3 &  r & w_1 & w_2 &   \\
    \hline
    f   &  0 &  0 &  0 &  0  &  0  &  0  &  0 & -1  & -1  & 0 \\
    s_1 &  1 &  1 &  1 &  1  &  0  &  0  &  0 &  0  &  0  & 2 \\
    w_1 &  1 &  0 & -1 &  0  &  0  &  0  &  0 &  1  &  0  & 1 \\
    w_2 &  0 &  1 &  0 &  0  &  0  &  0  & -1 &  0  &  1  & 1 \\
    s_2 &  1 &  0 &  0 &  0  &  1  &  0  &  0 &  0  &  0  & 3 \\
    s_3 & -1 &  1 &  0 &  0  &  0  &  1  &  0 &  0  &  0  & 1 \\
\end{array} 
\end{equation}
Den översta raden ska tolkas $f - w_1 - w_2 = 0$. Den är \emph{inte}
på kanonisk form, eftersom basvariablerna $w_1, w_2$ inte har
eliminerats. Men detta görs enkelt genom att addera deras rader till
målfunktionen
\begin{equation}
  \begin{array}{c|rrrrrrrrr|c}
        &  x & y & z  & s_1 & s_2 & s_3 &  r & w_1 & w_2 &   \\
        \hline
    f   &  1 & 1 & -1 &  0  &  0  &  0  & -1 &  0  &  0  & 2 \\
    s_1 &  1 & 1 &  1 &  1  &  0  &  0  &  0 &  0  &  0  & 2 \\
    w_1 &  1 & 0 & -1 &  0  &  0  &  0  &  0 &  1  &  0  & 1 \\
    w_2 &  0 & 1 &  0 &  0  &  0  &  0  & -1 &  0  &  1  & 1 \\
    s_2 &  1 & 0 &  0 &  0  &  1  &  0  &  0 &  0  &  0  & 3 \\
    s_3 & -1 & 1 &  0 &  0  &  0  &  1  &  0 &  0  &  0  & 1 \\
  \end{array}
\end{equation}
Möjliga inkommande basvariabler är $x$ och $y$. Vi väljer $x$. Eftersom
$b_1/A_{11} = 2, b_2/A_{21} = 1, b_4/A_{41} = 3$, så blir $w_1$ utgående
basvariabel. Den nya tablån blir
\begin{equation}
  \begin{array}{c|rrrrrrrrr|c}
        & x & y &  z & s_1 & s_2 & s_3 &  r & w_1 & w_2 &   \\
        \hline
    f   & 0 & 1 &  0 &  0  &  0  &  0  & -1 & -1  &  0  & 1 \\
    s_1 & 0 & 1 &  2 &  1  &  0  &  0  &  0 & -1  &  0  & 1 \\
    x   & 1 & 0 & -1 &  0  &  0  &  0  &  0 &  1  &  0  & 1 \\
    w_2 & 0 & 1 &  0 &  0  &  0  &  0  & -1 &  0  &  1  & 1 \\
    s_2 & 0 & 0 &  1 &  0  &  1  &  0  &  0 & -1  &  0  & 2 \\
    s_3 & 0 & 1 & -1 &  0  &  0  &  1  &  0 &  1  &  0  & 2 \\
  \end{array}
\end{equation}
Nu väljer vi $y$ som inkommande basvariabel, och möjliga utgående
basvariabler blir $s_1, w_2$. Eftersom målet är att få ner de
artificiella variablerna till 0 så borde man välja $w_2$ som utgående
variabel, men vi leker dum datanmaskin och väljer $s_1$ eftersom den är
först.
\begin{equation}
  \begin{array}{c|rrrrrrrrr|c}
        & x & y &  z & s_1 & s_2 & s_3 &  r & w_1 & w_2 &   \\
        \hline
    f   & 0 & 0 & -2 & -1  &  0  &  0  & -1 &  0  &  0  & 0 \\
    y   & 0 & 1 &  2 &  1  &  0  &  0  &  0 & -1  &  0  & 1 \\
    x   & 1 & 0 & -1 &  0  &  0  &  0  &  0 &  1  &  0  & 1 \\
    w_2 & 0 & 0 & -2 & -1  &  0  &  0  & -1 &  1  &  1  & 0 \\
    s_2 & 0 & 0 &  1 &  0  &  1  &  0  &  0 & -1  &  0  & 2 \\
    s_3 & 0 & 0 & -3 & -1  &  0  &  1  &  0 &  2  &  0  & 1 \\
  \end{array}
\end{equation}
Nu är alla $c_i \leq 0$, så vi har kommit fram till minimum. Vi ser,
från första raden, att $f = 0$ är minimum, och kan dra slutsatsen att
triangeln skär polyedern. Vi kan också läsa ut värdena på alla
variabler: $z = s_1 = r = w_1 = 0$ (icke-basvariabler), och $y = 1, x
= 1, w_2 = 0, s_2 = 2, s_3 = 1$. Man kan kontrollera att detta
verkligen är en lösning till de ursprungliga ekvationssystemet.

Slutsatsen är alltså att polyedern och triangeln skär varandra, och
att snittet innehåller punkten $(x, y, z) = (1, 1, 0)$.

\end{document}
