nx_cornell.nx source
↩ module page · 333 lines · 16757 B
1// nx_cornell.nx -- ★THE CORNELL BOX on the sovereign stack (field-evidence rung 3: the measured GI reference). A real
2// PATCH RADIOSITY solver -- Goral/Torrance/Greenberg/Battaile's own SIGGRAPH'84 method, for the very box they built:
3// red left wall, green right wall, white floor/ceiling/back, ceiling area light, two blocks. Patches gather B = E +
4// rho*Sum(F_ij*B_j) with point-to-patch form factors + BLOCK OCCLUSION (soft shadows); COLOR BLEEDING emerges from the
5// bounce, never painted. Axis-aligned block variant; coords shifted so the box centre is the origin (camera z=-camz).
6// All integer: F in Q16, radiosity B in milli, reflectance fx1024. license_tier: ORIGINAL
7import "nx_syscalls.nx"
8import "nx_trimesh.nx"
9import "nx_vecmath.nx"
10const COR_MAGIC_100000: i64 = 100000
11const COR_MAGIC_1024: i64 = 1024
12const COR_MAGIC_51472: i64 = 51472
13const COR_MAGIC_16384: i64 = 16384
14const COR_MAGIC_65536: i64 = 65536
15const COR_MAGIC_4096: i64 = 4096
16
17const COR_CAP: i64 = 700
18static COR_NP: i64
19static COR_CX: i64 // patch centers
20static COR_CY: i64
21static COR_CZ: i64
22static COR_NXA: i64 // unit normal components (-1/0/1)
23static COR_NYA: i64
24static COR_NZA: i64
25static COR_AR: i64 // patch area (units^2)
26static COR_RR: i64 // reflectance fx1024
27static COR_RG: i64
28static COR_RB: i64
29static COR_ER: i64 // emission (milli)
30static COR_EG: i64
31static COR_EB: i64
32static COR_BR: i64 // radiosity B (milli)
33static COR_BG2: i64
34static COR_BB: i64
35static COR_TR: i64 // scratch B'
36static COR_TG: i64
37static COR_TB: i64
38static COR_HUX: i64 // half-edge vectors (for the render quads)
39static COR_HUY: i64
40static COR_HUZ: i64
41static COR_HVX: i64
42static COR_HVY: i64
43static COR_HVZ: i64
44static COR_F: i64 // Q16 form factors, NP*NP
45
46func cor_isqrt(v: i64) -> i64 { return vm_isqrt(v) }
47func cor_alloc() -> i64 {
48 if COR_CX != 0 { return 0 }
49 COR_CX = sys_mmap(COR_CAP*8) as i64; COR_CY = sys_mmap(COR_CAP*8) as i64; COR_CZ = sys_mmap(COR_CAP*8) as i64
50 COR_NXA = sys_mmap(COR_CAP*8) as i64; COR_NYA = sys_mmap(COR_CAP*8) as i64; COR_NZA = sys_mmap(COR_CAP*8) as i64
51 COR_AR = sys_mmap(COR_CAP*8) as i64
52 COR_RR = sys_mmap(COR_CAP*8) as i64; COR_RG = sys_mmap(COR_CAP*8) as i64; COR_RB = sys_mmap(COR_CAP*8) as i64
53 COR_ER = sys_mmap(COR_CAP*8) as i64; COR_EG = sys_mmap(COR_CAP*8) as i64; COR_EB = sys_mmap(COR_CAP*8) as i64
54 COR_BR = sys_mmap(COR_CAP*8) as i64; COR_BG2 = sys_mmap(COR_CAP*8) as i64; COR_BB = sys_mmap(COR_CAP*8) as i64
55 COR_TR = sys_mmap(COR_CAP*8) as i64; COR_TG = sys_mmap(COR_CAP*8) as i64; COR_TB = sys_mmap(COR_CAP*8) as i64
56 COR_HUX = sys_mmap(COR_CAP*8) as i64; COR_HUY = sys_mmap(COR_CAP*8) as i64; COR_HUZ = sys_mmap(COR_CAP*8) as i64
57 COR_HVX = sys_mmap(COR_CAP*8) as i64; COR_HVY = sys_mmap(COR_CAP*8) as i64; COR_HVZ = sys_mmap(COR_CAP*8) as i64
58 COR_F = sys_mmap(COR_CAP*COR_CAP*8) as i64
59 return 0
60}
61// add a KxK grid of patches on the rectangle O + U*s + V*t. P = [x0,y0,z0, ux,uy,uz, vx,vy,vz, K, nxa,nya,nza,
62// rr,rg,rb, er,eg,eb] (bundled pointer -- dodges the 16-arg cap)
63func cor_addgrid(P: *i64) -> i64 {
64 let K: i64 = P[9]
65 let CXp: *i64 = COR_CX as *i64; let CYp: *i64 = COR_CY as *i64; let CZp: *i64 = COR_CZ as *i64
66 let NX: *i64 = COR_NXA as *i64; let NY: *i64 = COR_NYA as *i64; let NZ: *i64 = COR_NZA as *i64
67 let AR: *i64 = COR_AR as *i64
68 let RR: *i64 = COR_RR as *i64; let RG: *i64 = COR_RG as *i64; let RB: *i64 = COR_RB as *i64
69 let ER: *i64 = COR_ER as *i64; let EG: *i64 = COR_EG as *i64; let EB: *i64 = COR_EB as *i64
70 let HUX: *i64 = COR_HUX as *i64; let HUY: *i64 = COR_HUY as *i64; let HUZ: *i64 = COR_HUZ as *i64
71 let HVX: *i64 = COR_HVX as *i64; let HVY: *i64 = COR_HVY as *i64; let HVZ: *i64 = COR_HVZ as *i64
72 let ulen: i64 = cor_isqrt(P[3]*P[3] + P[4]*P[4] + P[5]*P[5])
73 let vlen: i64 = cor_isqrt(P[6]*P[6] + P[7]*P[7] + P[8]*P[8])
74 let parea: i64 = ulen*vlen/(K*K)
75 var i: i64 = 0
76 while i < K {
77 var j: i64 = 0
78 while j < K {
79 let m: i64 = COR_NP
80 CXp[m] = P[0] + P[3]*(2*i+1)/(2*K) + P[6]*(2*j+1)/(2*K)
81 CYp[m] = P[1] + P[4]*(2*i+1)/(2*K) + P[7]*(2*j+1)/(2*K)
82 CZp[m] = P[2] + P[5]*(2*i+1)/(2*K) + P[8]*(2*j+1)/(2*K)
83 NX[m] = P[10]; NY[m] = P[11]; NZ[m] = P[12]
84 AR[m] = parea
85 RR[m] = P[13]; RG[m] = P[14]; RB[m] = P[15]
86 ER[m] = P[16]; EG[m] = P[17]; EB[m] = P[18]
87 HUX[m] = P[3]/(2*K); HUY[m] = P[4]/(2*K); HUZ[m] = P[5]/(2*K)
88 HVX[m] = P[6]/(2*K); HVY[m] = P[7]/(2*K); HVZ[m] = P[8]/(2*K)
89 COR_NP = m + 1
90 j = j + 1
91 }
92 i = i + 1
93 }
94 return 0
95}
96// ★build the canonical scene (shifted: x' = x-278 [-276..276], y' = y-273 [-273..275], z 0..559; camera at z=-camz).
97// Layout (for the gate): floor 0..99 (K=10, i over x, j over z) | ceiling 100..199 | BACK 200..343 (K=12, i over x,
98// j over y) | RED left 344..443 | GREEN right 444..543 | light 544..559 | tall block 560..604 | short block 605..649.
99func cornell_build() -> i64 {
100 cor_alloc()
101 cor_boxes_init()
102 COR_NP = 0
103 let P: *i64 = sys_mmap(19*8) as *i64
104 let WH: i64 = 740 // neutral white reflectance (EXACTLY r=g=b -> direct r-g == 0)
105 // floor y=-273
106 P[0]=0-276; P[1]=0-273; P[2]=0; P[3]=552; P[4]=0; P[5]=0; P[6]=0; P[7]=0; P[8]=559
107 P[9]=10; P[10]=0; P[11]=1; P[12]=0; P[13]=WH; P[14]=WH; P[15]=WH; P[16]=0; P[17]=0; P[18]=0
108 cor_addgrid(P)
109 // ceiling y=275
110 P[1]=275; P[10]=0; P[11]=0-1; P[12]=0
111 cor_addgrid(P)
112 // back wall z=559 (K=12; i over x, j over y) -- the bleeding is measured HERE
113 P[0]=0-276; P[1]=0-273; P[2]=559; P[3]=552; P[4]=0; P[5]=0; P[6]=0; P[7]=548; P[8]=0
114 P[9]=12; P[10]=0; P[11]=0; P[12]=0-1
115 cor_addgrid(P)
116 // RED left wall x=-276 (canonical red ~ (0.63,0.06,0.05))
117 P[0]=0-276; P[1]=0-273; P[2]=0; P[3]=0; P[4]=0; P[5]=559; P[6]=0; P[7]=548; P[8]=0
118 P[9]=10; P[10]=1; P[11]=0; P[12]=0
119 P[13]=645; P[14]=62; P[15]=52
120 cor_addgrid(P)
121 // GREEN right wall x=276 (canonical green ~ (0.12,0.45,0.10))
122 P[0]=276; P[10]=0-1
123 P[13]=120; P[14]=460; P[15]=104
124 cor_addgrid(P)
125 // light: ceiling area lamp (130x105) 1 unit below the ceiling, white emitter
126 P[0]=0-65; P[1]=274; P[2]=227; P[3]=130; P[4]=0; P[5]=0; P[6]=0; P[7]=0; P[8]=105
127 P[9]=4; P[10]=0; P[11]=0-1; P[12]=0
128 P[13]=740; P[14]=740; P[15]=740
129 P[16]=COR_MAGIC_100000; P[17]=COR_MAGIC_100000; P[18]=COR_MAGIC_100000
130 cor_addgrid(P)
131 P[16]=0; P[17]=0; P[18]=0
132 P[13]=WH; P[14]=WH; P[15]=WH
133 // tall block (left-back): x -180..-60, y -273..57, z 250..400 -- 5 faces, K=3
134 P[9]=3
135 P[0]=0-180; P[1]=0-273; P[2]=250; P[3]=120; P[4]=0; P[5]=0; P[6]=0; P[7]=330; P[8]=0; P[10]=0; P[11]=0; P[12]=0-1
136 cor_addgrid(P) // front z=250
137 P[2]=400; P[12]=1
138 cor_addgrid(P) // back z=400
139 P[0]=0-180; P[2]=250; P[3]=0; P[4]=0; P[5]=150; P[6]=0; P[7]=330; P[8]=0; P[10]=0-1; P[12]=0
140 cor_addgrid(P) // left x=-180
141 P[0]=0-60; P[10]=1
142 cor_addgrid(P) // right x=-60
143 P[0]=0-180; P[1]=57; P[2]=250; P[3]=120; P[4]=0; P[5]=0; P[6]=0; P[7]=0; P[8]=150; P[10]=0; P[11]=1
144 cor_addgrid(P) // top y=57
145 // short block (right-front): x 30..170, y -273..-108, z 80..230 -- 5 faces, K=3
146 P[0]=30; P[1]=0-273; P[2]=80; P[3]=140; P[4]=0; P[5]=0; P[6]=0; P[7]=165; P[8]=0; P[10]=0; P[11]=0; P[12]=0-1
147 cor_addgrid(P) // front z=80
148 P[2]=230; P[12]=1
149 cor_addgrid(P) // back z=230
150 P[0]=30; P[2]=80; P[3]=0; P[4]=0; P[5]=150; P[6]=0; P[7]=165; P[8]=0; P[10]=0-1; P[12]=0
151 cor_addgrid(P) // left x=30
152 P[0]=170; P[10]=1
153 cor_addgrid(P) // right x=170
154 P[0]=30; P[1]=0-108; P[2]=80; P[3]=140; P[4]=0; P[5]=0; P[6]=0; P[7]=0; P[8]=150; P[10]=0; P[11]=1
155 cor_addgrid(P) // top y=-108
156 return COR_NP
157}
158// segment (x0..x1) vs one inset AABB (slab, t in ~[2%,98%] of the segment). 1 = blocked.
159func cor_seg_box(x0: i64, y0: i64, z0: i64, x1: i64, y1: i64, z1: i64, B: *i64) -> i64 {
160 var tmin: i64 = 20 // epsilon: skip endpoints (patches lie ON the boxes)
161 var tmax: i64 = 1004
162 var axis: i64 = 0
163 while axis < 3 {
164 var o: i64 = x0; var d: i64 = x1-x0; var lo: i64 = B[0]; var hi: i64 = B[3]
165 if axis == 1 { o = y0; d = y1-y0; lo = B[1]; hi = B[4] }
166 if axis == 2 { o = z0; d = z1-z0; lo = B[2]; hi = B[5] }
167 if d == 0 {
168 if o < lo { return 0 }
169 if o > hi { return 0 }
170 } else {
171 var t1: i64 = (lo-o)*COR_MAGIC_1024/d
172 var t2: i64 = (hi-o)*COR_MAGIC_1024/d
173 if t1 > t2 { let tt: i64 = t1; t1 = t2; t2 = tt }
174 if t1 > tmin { tmin = t1 }
175 if t2 < tmax { tmax = t2 }
176 }
177 axis = axis + 1
178 }
179 if tmax >= tmin { return 1 }
180 return 0
181}
182// occluder boxes as STATE (dynamic geometry: cornell_move_tall shifts box 1 with its patches)
183static COR_BOX1: i64
184static COR_BOX2: i64
185func cor_boxes_init() -> i64 {
186 if COR_BOX1 == 0 { COR_BOX1 = sys_mmap(48) as i64; COR_BOX2 = sys_mmap(48) as i64 }
187 let B1: *i64 = COR_BOX1 as *i64
188 B1[0]=0-179; B1[1]=0-272; B1[2]=251; B1[3]=0-61; B1[4]=56; B1[5]=399
189 let B2: *i64 = COR_BOX2 as *i64
190 B2[0]=31; B2[1]=0-272; B2[2]=81; B2[3]=169; B2[4]=0-109; B2[5]=229
191 return 0
192}
193// is the segment i->j occluded by either block? (blocks inset by 1 so a face's own patches don't self-occlude)
194func cor_occluded(x0: i64, y0: i64, z0: i64, x1: i64, y1: i64, z1: i64) -> i64 {
195 if cor_seg_box(x0,y0,z0,x1,y1,z1, COR_BOX1 as *i64) == 1 { return 1 }
196 if cor_seg_box(x0,y0,z0,x1,y1,z1, COR_BOX2 as *i64) == 1 { return 1 }
197 return 0
198}
199// ---- DYNAMIC SCENE EVENTS (R3): move the light / move the tall block ----
200// layout constants: light patches 544..559, tall-block patches 560..604 (see cornell_build comment)
201func cornell_move_light(dx: i64, dz: i64) -> i64 {
202 let CXp: *i64 = COR_CX as *i64
203 let CZp: *i64 = COR_CZ as *i64
204 var i: i64 = 544
205 while i < 560 { CXp[i] = CXp[i] + dx; CZp[i] = CZp[i] + dz; i = i + 1 }
206 return 0
207}
208func cornell_move_tall(dx: i64) -> i64 {
209 let CXp: *i64 = COR_CX as *i64
210 var i: i64 = 560
211 while i < 605 { CXp[i] = CXp[i] + dx; i = i + 1 }
212 let B1: *i64 = COR_BOX1 as *i64
213 B1[0] = B1[0] + dx
214 B1[3] = B1[3] + dx
215 return 0
216}
217// one Q16 point-to-patch form factor with occlusion: F_ij = cosi*cosj*Aj / (pi*(r2+Aj)), 0 if occluded.
218func cor_ff_pair(i: i64, j: i64) -> i64 {
219 if j == i { return 0 }
220 let CXp: *i64 = COR_CX as *i64; let CYp: *i64 = COR_CY as *i64; let CZp: *i64 = COR_CZ as *i64
221 let NX: *i64 = COR_NXA as *i64; let NY: *i64 = COR_NYA as *i64; let NZ: *i64 = COR_NZA as *i64
222 let AR: *i64 = COR_AR as *i64
223 let dx: i64 = CXp[j]-CXp[i]; let dy: i64 = CYp[j]-CYp[i]; let dz: i64 = CZp[j]-CZp[i]
224 let doti: i64 = NX[i]*dx + NY[i]*dy + NZ[i]*dz // >0 if j in front of i
225 let dotj: i64 = 0 - (NX[j]*dx + NY[j]*dy + NZ[j]*dz) // >0 if i in front of j
226 if doti <= 0 { return 0 }
227 if dotj <= 0 { return 0 }
228 let r2: i64 = dx*dx + dy*dy + dz*dz
229 let len: i64 = cor_isqrt(r2) + 1
230 let ci: i64 = doti*COR_MAGIC_1024/len
231 let cj: i64 = dotj*COR_MAGIC_1024/len
232 if cor_occluded(CXp[i],CYp[i],CZp[i], CXp[j],CYp[j],CZp[j]) == 1 { return 0 }
233 let t: i64 = ci*cj/COR_MAGIC_1024
234 let u: i64 = t*AR[j]
235 let den: i64 = COR_MAGIC_51472*((r2 + AR[j])/COR_MAGIC_1024) + 51 // 16*pi*(r2+Aj) in /COR_MAGIC_1024 scale
236 var f: i64 = u*COR_MAGIC_1024/den
237 if f > COR_MAGIC_16384 { f = COR_MAGIC_16384 } // clamp F <= 0.25
238 return f
239}
240// ★precompute ALL Q16 form factors (the cached "surface" of the scene). Returns pairs computed.
241func cornell_formfactors() -> i64 {
242 let F: *i64 = COR_F as *i64
243 let NP: i64 = COR_NP
244 var i: i64 = 0
245 while i < NP {
246 var j: i64 = 0
247 while j < NP { F[i*NP+j] = cor_ff_pair(i, j); j = j + 1 }
248 i = i + 1
249 }
250 return NP*NP
251}
252// ★INCREMENTAL invalidation (the Surface-Cache move): recompute ONLY pairs touching patches [from,to).
253// Correct when the moved patches' geometry changed but occluders did not. Returns pairs recomputed.
254func cornell_ff_partial(from: i64, to: i64) -> i64 {
255 let F: *i64 = COR_F as *i64
256 let NP: i64 = COR_NP
257 var work: i64 = 0
258 var i: i64 = from
259 while i < to {
260 var j: i64 = 0
261 while j < NP {
262 F[i*NP+j] = cor_ff_pair(i, j)
263 F[j*NP+i] = cor_ff_pair(j, i)
264 work = work + 2
265 j = j + 1
266 }
267 i = i + 1
268 }
269 return work
270}
271// ★solve radiosity: B = E; `iters` gather sweeps of B' = E + rho*Sum(F*B)/Q16. iters=1 -> DIRECT only; 3 -> +2 bounces.
272func cornell_solve(iters: i64) -> i64 {
273 let ER: *i64 = COR_ER as *i64; let EG: *i64 = COR_EG as *i64; let EB: *i64 = COR_EB as *i64
274 let RR: *i64 = COR_RR as *i64; let RG: *i64 = COR_RG as *i64; let RB: *i64 = COR_RB as *i64
275 let BR: *i64 = COR_BR as *i64; let BG: *i64 = COR_BG2 as *i64; let BB: *i64 = COR_BB as *i64
276 let TR: *i64 = COR_TR as *i64; let TG: *i64 = COR_TG as *i64; let TB: *i64 = COR_TB as *i64
277 let F: *i64 = COR_F as *i64
278 let NP: i64 = COR_NP
279 var i: i64 = 0
280 while i < NP { BR[i]=ER[i]; BG[i]=EG[i]; BB[i]=EB[i]; i=i+1 }
281 var it: i64 = 0
282 while it < iters {
283 i = 0
284 while i < NP {
285 var gr: i64 = 0
286 var gg: i64 = 0
287 var gb: i64 = 0
288 var j: i64 = 0
289 while j < NP {
290 let f: i64 = F[i*NP+j]
291 if f != 0 { gr = gr + f*BR[j]; gg = gg + f*BG[j]; gb = gb + f*BB[j] }
292 j = j + 1
293 }
294 TR[i] = ER[i] + RR[i]*(gr/COR_MAGIC_65536)/COR_MAGIC_1024
295 TG[i] = EG[i] + RG[i]*(gg/COR_MAGIC_65536)/COR_MAGIC_1024
296 TB[i] = EB[i] + RB[i]*(gb/COR_MAGIC_65536)/COR_MAGIC_1024
297 i = i + 1
298 }
299 i = 0
300 while i < NP { BR[i]=TR[i]; BG[i]=TG[i]; BB[i]=TB[i]; i=i+1 }
301 it = it + 1
302 }
303 return 0
304}
305func cor_np() -> i64 { return COR_NP }
306func cor_b(i: i64, out: *i64) -> i64 { let BR: *i64=COR_BR as *i64; let BG: *i64=COR_BG2 as *i64; let BB: *i64=COR_BB as *i64; out[0]=BR[i]; out[1]=BG[i]; out[2]=BB[i]; return 0 }
307func cor_center(i: i64, out: *i64) -> i64 { let CXp: *i64=COR_CX as *i64; let CYp: *i64=COR_CY as *i64; let CZp: *i64=COR_CZ as *i64; out[0]=CXp[i]; out[1]=CYp[i]; out[2]=CZp[i]; return 0 }
308// emit the patches as UNLIT vertex-coloured quads (render with smooth==3). gain: display = B*255/gain clamp.
309func cornell_tomesh(gain: i64) -> i64 {
310 let CXp: *i64 = COR_CX as *i64; let CYp: *i64 = COR_CY as *i64; let CZp: *i64 = COR_CZ as *i64
311 let NX: *i64 = COR_NXA as *i64; let NY: *i64 = COR_NYA as *i64; let NZ: *i64 = COR_NZA as *i64
312 let BR: *i64 = COR_BR as *i64; let BG: *i64 = COR_BG2 as *i64; let BB: *i64 = COR_BB as *i64
313 let HUX: *i64 = COR_HUX as *i64; let HUY: *i64 = COR_HUY as *i64; let HUZ: *i64 = COR_HUZ as *i64
314 let HVX: *i64 = COR_HVX as *i64; let HVY: *i64 = COR_HVY as *i64; let HVZ: *i64 = COR_HVZ as *i64
315 tm_reset()
316 var i: i64 = 0
317 while i < COR_NP {
318 var r: i64 = BR[i]*255/gain
319 var g: i64 = BG[i]*255/gain
320 var b: i64 = BB[i]*255/gain
321 if r > 255 { r = 255 } if g > 255 { g = 255 } if b > 255 { b = 255 }
322 let col: i64 = r + g*256 + b*COR_MAGIC_65536
323 let nx: i64 = NX[i]*COR_MAGIC_4096; let ny: i64 = NY[i]*COR_MAGIC_4096; let nz: i64 = NZ[i]*COR_MAGIC_4096
324 let v0: i64 = tm_vert_n(CXp[i]-HUX[i]-HVX[i], CYp[i]-HUY[i]-HVY[i], CZp[i]-HUZ[i]-HVZ[i], nx,ny,nz)
325 let v1: i64 = tm_vert_n(CXp[i]+HUX[i]-HVX[i], CYp[i]+HUY[i]-HVY[i], CZp[i]+HUZ[i]-HVZ[i], nx,ny,nz)
326 let v2: i64 = tm_vert_n(CXp[i]+HUX[i]+HVX[i], CYp[i]+HUY[i]+HVY[i], CZp[i]+HUZ[i]+HVZ[i], nx,ny,nz)
327 let v3: i64 = tm_vert_n(CXp[i]-HUX[i]+HVX[i], CYp[i]-HUY[i]+HVY[i], CZp[i]-HUZ[i]+HVZ[i], nx,ny,nz)
328 tm_vcol(v0,col); tm_vcol(v1,col); tm_vcol(v2,col); tm_vcol(v3,col)
329 tm_quad(v0,v1,v2,v3,col)
330 i = i + 1
331 }
332 return 0
333}