Raytracing med simplexalgoritmen

Niels Möller <nisse@lysator.liu.se> 4 juni 1996

Innehåll

Bakgrund 2 Simplexanpassad problemformulering 2 Simplexalgoritmen 4 Exempel 5

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 \overline{n} och sitt avstånd c till origo:

P = \{\overline{r} = (x, y, z) : \overline{r} · \overline{n} = c\}  (1)

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

for t in triangles:
   for p in polyeders:
      if intersectp(p, t):
         p1, p2 = split(p, t)
         polyeders.replace(p, p1, p2)

Här betyder split att låta det plan som triangeln ligger i skära polyderern i två delar, vilket är rätt enkelt. Det kluriga är 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.

Simplexanpassad problemformulering

Vi vill undersöka snittet mellan en triangel och en konvex polyeder. Snittet ges av lösningsmängden till

\overline{r} · \overline{n} = c  (2)
\overline{r} · \overline{a_{i}} ≤ α_{i} , för 1 ≤ i ≤ k  (3)
\overline{r} · \overline{b_{i}} ≥ β_{i} , för 1 ≤ i ≤ l  (4)

Här ska samtliga c, α_{i}, β_{i} ≥ 0. Ekvation (2) 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 slackvariabler, för att kunna ersätta olikheterna med likheter. Då ersätts ekvationerna (3) och (4) av

\overline{r} · \overline{a_{i}} + s_{i} = α_{i} , för 1 ≤ i ≤ k  (5)
\overline{r} · \overline{b_{i}} − r_{i} = β_{i} , för 1 ≤ i ≤ l  (6)

Även r_{i}, s_{i} ≥ 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 tillåtna baslösningar, vilket betyder att man sätter alla varibler utom m stycken till 0 (så kallade icke-basvariabler) och löser ut de övriga m 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} = α_{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 artificiella variabler i ekvationerna (2) och (6). Då får vi ett nytt ekvationssystem,

\overline{r} · \overline{n} + w_{0} = c  (7)
\overline{r} · \overline{a_{i}} + s_{i} = α_{i}, 1 ≤ i ≤ k  (8)
\overline{r} · \overline{b_{i}} − r_{i} + w_{i} = β_{i}, 1 ≤ i ≤ l  (9)
x, y, z, s_{i}, r_{i}, w_{i} ≥ 0  (10)

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

f = ∑_{0}^{l} w_{i}  (11)

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} = α_{i} för 1 ≤ i ≤ k, w_{0} = c, w_{i} = β_{i} för 1 ≤ i ≤ 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 (11), och flyttar över konstanterna, får vi

f + \overline{r} · (\overline{n} + ∑_{1}^{l} \overline{b_{i}}) − ∑_{1}^{l} r_{i} = c + ∑_{1}^{l} β_{i}  (12)

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 (12), så tar vi någon av de variabler vars koefficient är positiv, och ökar den. Detta är vår 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

• 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 < ϵ.

• 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 ∑ w_{i} kan aldrig bli negativ. Men konstiga avrundningsfel skulle möjligen kunna stoppa simplexalgoritmen här i alla fall.

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 × 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

min f då f + c^{T} x = b_{0}  (13)
Ax = b  (14)
x_{i} ≥ 0 , för 1 ≤ i ≤ n  (15)

När vi arbetar med problemet vill vi hela tiden ha det på kanonisk form, vilket betyder att basvariablerna ska vara eliminerade ur (13), och att motsvarande kolonner i A ska bilda en enhetsmatris. Med andra ord ska A_{kB_{j}} = δ_{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 exakt 0 och 1.

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

• Välj ny inkommande basvariabel.

• Välj utgående basvariabel.

• Gör radoperationer i A, b och c tills problemet är på kanonisk form (för den nya uppsättningen basvariabler).

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

1. Om alla c_{i} ≤ 0, stoppa.

2. Välj ett r, med c_{r} > 0 (med eller utan Blands regel).

3. Om alla A_{jr} ≤ 0, stoppa.

4. Tag minimum av {b_{j}/A_{jr}|A_{jr} > 0}, och låt k vara det lägsta index som ger minimum.

5. Sätt B′_{j} = B_{j} för j ≠ k, B′_{k} = r.

6. Sätt A′_{k} = A_{k}/A_{kr} och b′_{k} = b_{k}/A_{kr}. (A_{k} betecknar här rad k i A).

7. För j ≠ k, sätt A′_{j} = A_{j} − A_{jr}A_{k} och b′_{j} = b_{j} − A_{jr}b_{k}.

8. Sätt c′ = c − c_{r}A_{k} och b′_{0} = b_{0} − c_{r}b_{k}

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

Exempel

Vi vill undersöka om polyedern

x + y + z ≤ 2  (16)
x, y, z ≥ 0  (17)

skär triangeln

x − z = 1  (18)
y ≥ 1  (19)
x ≤ 3  (20)
−x + y ≤ 1  (21)

Inför vi slackvariabler i olikheterna får vi ekvationssystemet

x + y + z + s_{1} = 2  (22)
x − z = 1  (23)
y − r = 1  (24)
x + s_{2} = 3  (25)
−x + y + s_{3} = 1  (26)
x, y, z, s_{1}, s_{2}, s_{3}, r ≥ 0  (27)

Vi saknar basvariabler i ekvationerna (23) och (24), och behöver därför införa två artificiella variabler w_{1}, w_{2}.

\begin{aligned}& x & y & z & s_{1} & s_{2} & s_{3} & r & w_{1} & w_{2} \\ 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{aligned}  (28)

Den översta raden ska tolkas f − w_{1} − w_{2} = 0. Den är 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{aligned}& x & y & z & s_{1} & s_{2} & s_{3} & r & w_{1} & w_{2} \\ 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{aligned}  (29)

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{aligned}& x & y & z & s_{1} & s_{2} & s_{3} & r & w_{1} & w_{2} \\ 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{aligned}  (30)

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{aligned}& x & y & z & s_{1} & s_{2} & s_{3} & r & w_{1} & w_{2} \\ 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{aligned}  (31)

Nu är alla c_{i} ≤ 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 (ickebasvariabler), 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).