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