#include <stdio.h>
#include <stdlib.h>
+#include <math.h>
#include <SDL/SDL.h>
#include "cgmath/cgmath.h"
static void cleanup(void);
static void update_proj_vp(void);
static void draw(void);
-static int quadric(float *obj, int x, int y, float *retz);
+static int quadric(float *obj, cgm_vec4 *p);
+static int calc_w(float *obj, cgm_vec4 *p);
static void proc_event(SDL_Event *ev);
static int quit, resize_pending;
static float win_aspect;
static SDL_Surface *fbsurf;
static int frame;
+static int mouse_x, mouse_y;
+static int bnstate[8];
+static unsigned long start_msec;
#define SDLFLAGS SDL_SWSURFACE | SDL_RESIZABLE
+static float cam_theta, cam_phi = 0.5, cam_dist = 4;
static float modelview[16], projection[16], viewport[16];
+static float invmodelview[16], invproj[16];
+static float mvadj[16], mvadjt[16];
static float mvpmat[16], mvpadj[16], mvpadjt[16];
static float objmat[16];
SDL_Quit();
return 1;
}
+ SDL_WM_SetCaption("algebraic surface rasterization", 0);
win_width = 800;
win_height = 600;
win_aspect = 1.33333333333;
return 1;
}
+ start_msec = SDL_GetTicks();
for(;;) {
while(SDL_PollEvent(&ev)) {
proc_event(&ev);
return 0;
}
+static const float sphere[] = {
+ 1, 0, 0, 0,
+ 0, 1, 0, 0,
+ 0, 0, 1, 0,
+ 0, 0, 0, -1
+};
+static const float cone[] = {
+ 1, 0, 0, 0,
+ 0, 1, 0, 0,
+ 0, 0, -1, 0,
+ 0, 0, 0, 0
+};
+
+
static int init(void)
{
- cgm_mzero(objmat);
- objmat[0] = objmat[5] = objmat[10] = 1;
- objmat[15] = -1;
-
+ memcpy(objmat, sphere, sizeof objmat);
return 0;
}
static void update_proj_vp(void)
{
cgm_mperspective(projection, cgm_deg_to_rad(50), win_aspect, 0.5, 500.0f);
+ cgm_mcopy(invproj, projection);
+ cgm_minverse(invproj);
cgm_mtranslation(viewport, 1, 1, 0);
cgm_mscale(viewport, win_width * 0.5f, win_height * 0.5f, 1);
{
int i, j;
uint32_t *fbptr;
- float objp[16];
+ float objp[16], objv[16];
+ static const cgm_vec3 lpos = {-6, 10, 10};
+ cgm_vec3 ldir, hvec;
+ float t;
+ unsigned long msec;
+
+ msec = SDL_GetTicks() - start_msec;
+ t = (float)msec / 1000.0f;
if(SDL_MUSTLOCK(fbsurf)) {
SDL_LockSurface(fbsurf);
}
- fbptr = fbsurf->pixels;
- cgm_mtranslation(modelview, 0, 0, 5);
+ cgm_midentity(modelview);
+ /* view matrix */
+ cgm_mpretranslate(modelview, 0, -0.5, -cam_dist);
+ cgm_mprerotate_x(modelview, cam_phi);
+ cgm_mprerotate_y(modelview, cam_theta);
+
+ ldir = lpos; cgm_vmul_m4v3(&ldir, modelview);
+ cgm_vnormalize(&ldir);
+
+
+ /* model matrix */
+ cgm_mpretranslate(modelview, 0, fabs(cos(t)), 0);
+
+ cgm_mcopy(invmodelview, modelview);
+ cgm_minverse(invmodelview);
cgm_mcopy(mvpmat, modelview);
cgm_mmul(mvpmat, projection);
cgm_mmul(mvpmat, viewport);
/* calculate the adjoint, and the adjoint-transpose */
- cgm_mcopy(mvpadj, mvpmat); cgm_minverse(mvpadj);
+ cgm_mcopy(mvpadj, mvpmat); cgm_mcofmatrix(mvpadj); cgm_mtranspose(mvpadj);
cgm_mcopy(mvpadjt, mvpadj); cgm_mtranspose(mvpadjt);
- /* transform object */
+ cgm_mcopy(mvadj, modelview); cgm_mcofmatrix(mvadj); cgm_mtranspose(mvadj);
+ cgm_mcopy(mvadjt, mvadj); cgm_mtranspose(mvadjt);
+
+ /* transform object to screen space */
cgm_mcopy(objp, mvpadj);
cgm_mmul(objp, objmat);
cgm_mmul(objp, mvpadjt);
+ /* transform object to view space */
+ cgm_mcopy(objv, mvadj);
+ cgm_mmul(objv, objmat);
+ cgm_mmul(objv, mvadjt);
+
+ hvec.x = ldir.x * 0.5f;
+ hvec.y = ldir.y * 0.5f;
+ hvec.z = (ldir.z + 1.0f) * 0.5f;
+ cgm_vnormalize(&hvec);
+
+ fbptr = (uint32_t*)fbsurf->pixels + (win_height - 1) * win_width;
for(i=0; i<win_height; i++) {
for(j=0; j<win_width; j++) {
- float z;
- if(quadric(objp, j, i, &z)) {
- *fbptr++ = 0xff0000;
+ cgm_vec4 p;
+ cgm_vec3 n;
+ float ndotl, ndoth, spec;
+ int r, g, b;
+ float u, v;
+
+ p.x = j;
+ p.y = i;
+ if(quadric(objp, &p)) {
+ /* transform back to local space */
+ p.x = p.x * 2.0f / win_width - 1.0f;
+ p.y = p.y * 2.0f / win_height - 1.0f;
+ p.w = 1.0f;
+ cgm_wmul_m4v4(&p, invproj);
+
+ /* compute the normal in view space */
+ n.x = p.x * objv[0] + p.y * objv[1] + p.z * objv[2] + p.w * objv[3];
+ n.y = p.x * objv[4] + p.y * objv[5] + p.z * objv[6] + p.w * objv[7];
+ n.z = p.x * objv[8] + p.y * objv[9] + p.z * objv[10] + p.w * objv[11];
+ cgm_vnormalize(&n);
+
+ /* transofrm p to local space to compute uv mapping */
+ cgm_wmul_m4v4(&p, invmodelview);
+ cgm_wscale(&p, 1.0f / p.w);
+ u = atan2(p.z, p.x) * 0.5f / CGM_PI + CGM_PI/2;
+ v = acos(p.y) / CGM_PI;
+
+ ndotl = cgm_vdot(&n, &ldir);
+ if(ndotl < 0) ndotl = 0;
+ ndoth = cgm_vdot(&n, &hvec);
+ if(ndoth < 0) ndoth = 0;
+ spec = pow(ndoth, 60.0f);
+
+ u *= 16.0f;
+ v *= 8.0f;
+ u = u - floor(u);
+ v = v - floor(v);
+ if(u < 0.15 || v < 0.15) {
+ r = 64;
+ g = 128;
+ b = 255;
+ } else {
+ r = 255;
+ g = 80;
+ b = 80;
+ }
+
+ r = (int)(ndotl * r + spec * 192.0f) + r / 16;
+ g = (int)(ndotl * g + spec * 192.0f) + g / 16;
+ b = (int)(ndotl * b + spec * 192.0f) + b / 16;
+ if(r > 255) r = 255;
+ if(g > 255) g = 255;
+ if(b > 255) b = 255;
+
+ fbptr[j] = (r << 16) | (g << 8) | b;
} else {
- *fbptr++ = 0;
+ fbptr[j] = 0;
}
}
+ fbptr -= win_width;
}
if(SDL_MUSTLOCK(fbsurf)) {
SDL_Flip(fbsurf);
}
-static int quadric(float *obj, int xpix, int ypix, float *retz)
+static int quadric(float *obj, cgm_vec4 *p)
{
float x, y, a, b, c, d, sqrtd, a2, z0, z1;
- x = (float)xpix;
- y = (float)ypix;
+ x = p->x;
+ y = p->y;
a = obj[10];
b = (x * obj[2] + y * obj[6] + obj[11]) * 2.0f;
c = x * x * obj[0] + y * y * obj[5] + (x * y * obj[1] + x * obj[3] +
y * obj[7]) * 2.0f + obj[15];
- d = b * b - 4 * a * c;
+ d = b * b - 4.0f * a * c;
if(d <= 0.0f) return 0;
sqrtd = sqrt(d);
if(z0 < 0.0f) z0 = z1;
if(z1 < 0.0f) z1 = z0;
- *retz = z0 < z1 ? z0 : z1;
+ p->z = z0 < z1 ? z0 : z1;
+ return 1;
+}
+
+static int calc_w(float *obj, cgm_vec4 *p)
+{
+ float x, y, z, a, b, c, d, sqrtd, a2, w0, w1;
+
+ x = p->x;
+ y = p->y;
+ z = p->z;
+
+ a = obj[15];
+ b = (x * obj[3] + y * obj[7] + z * obj[11]) * 2.0f;
+ c = x * x * obj[0] + y * y * obj[5] + z * z * obj[10] +
+ (x * y * obj[1] + x * z * obj[2] + y * z * obj[6]) * 2.0f;
+
+ d = b * b - 4.0f * a * c;
+ if(d <= 0.0f) return 0;
+
+ sqrtd = sqrt(d);
+ a2 = a * 2.0f;
+
+ w0 = (-b + sqrtd) / a2;
+ w1 = (-b - sqrtd) / a2;
+
+ p->w = w0 < w1 ? w0 : w1; /* ???? */
return 1;
}
static void proc_event(SDL_Event *ev)
{
+ int bidx, dx, dy;
+
switch(ev->type) {
case SDL_KEYDOWN:
if(ev->key.keysym.sym == SDLK_ESCAPE) {
resize_pending = 1;
}
break;
+
+ case SDL_MOUSEBUTTONDOWN:
+ case SDL_MOUSEBUTTONUP:
+ bidx = ev->button.button - SDL_BUTTON_LEFT;
+ bnstate[bidx] = ev->button.state == SDL_PRESSED ? 1 : 0;
+ mouse_x = ev->button.x;
+ mouse_y = ev->button.y;
+ break;
+
+ case SDL_MOUSEMOTION:
+ dx = ev->motion.xrel;
+ dy = ev->motion.yrel;
+
+ if((dx | dy)) {
+ if(bnstate[0]) {
+ cam_theta += dx * 0.01;
+ cam_phi += dy * 0.01;
+ if(cam_phi < -CGM_PI/2) cam_phi = -CGM_PI/2;
+ if(cam_phi > CGM_PI/2) cam_phi = CGM_PI/2;
+ }
+ if(bnstate[2]) {
+ cam_dist += dy * 0.1;
+ if(cam_dist < 0) cam_dist = 0;
+ }
+ }
+ break;
}
}
projects / algsurf / commitdiff