/* Assignment 6 - wavelet reconstruction of landscape height field. * Its derived from the fractal assignment but now uses the 'h' and 'l' keys * to show Higher & Lower levels of detail from the wavelet reconstruction. * The goal is to get a general understanding of hierarchical wavelet ideas. * Since the wavelet operations work directly on the simg array we recopy * from ht_bytes on every lev_detail change and recompute both forward and * reverse transforms. Your knockout of wavelet terms will be in between * the forward and reverse transforms. Although simg is allocated as 1D * storage you will need to access as 2D like simg[y*ht_x+x]. You only * need to turn in the knockout code as 1) below but please send in any * other interesting experiments as well. * * Your job: * is now VERY easy and is experimentation and observation than real work. * 1) XXX implement the knockout of coefficients according to lev_detail. * That is, if lev_detal == 0 then everything is knocked out. * 2) XXX observe the effects of wavelet prediction define PRED or NOPRED * 3) XXX try differentially zapping vert, horizonal or diag wavelet terms * or taking out low rather than high coefficients * * To compile on linux: * gcc -Wall asn05.c -lglut -lGLU -lGL -lXmu -lXi -lX11 -lm -L/usr/X11R6/lib * * To compile on windows add: * glut32.lib glu32.lib opengl32.lib * * To Compile on MacOS X add: * -framework GLUT -framework OpenGL * -framework Foundation -framework CoreServices * -framework ApplicationServices -framework Carbon * * Hit the ESC (escape) key to exit */ /* required headers */ #include #include #include #include #include #include /* Constants */ #define WINDOW_X 1024 #define WINDOW_Y 768 #define PI 3.14159265358979323846 #define HF_ROW 1024 #define HF_COL 1024 /* containers */ typedef struct height_field_point { float height ; float color[3] ; float normal[3] ; } height_field_point ; /* Global state variables */ int mouse_button, mouse_state; int viewport_x, viewport_y ; int go = 1 ; // for the timer callback float cam_r = 500, cam_theta = PI, cam_phi = -PI/4.0 ; // radius, longitude, lat float cam_dtheta = 0, cam_dphi = 0, cam_dr = 0 ; float cam_d = PI/64.0 ; int ticks = 0 ; char draw_axes_flag = 0 ; int ht_x, ht_y; // dimension of ht map unsigned char *ht_bytes; // the actual ht map byte data short *simg; // short array used for wavelet computations int lev_detail; /* height field array, note the extra row/col for easy division */ height_field_point hf[HF_ROW+1][HF_COL+1] ; char mode = 0 ; char polymode = 0 ; /* elementary 3 component vector operations ~by order of class presentation */ /* return the length (ie. magnitude) of a 3 component isotropic vector */ float length(float a[3]) { return( sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]) ) ; } /* modify the vector a to normalize its length to 1.0 */ void normalize(float a[3]) { float len = length(a) ; if( len == 0 ) { fprintf(stderr,"attempt to divide by 0 in normalize()\n") ; return ; } a[0] /= len ; a[1] /= len ; a[2] /= len ; } /* add two vectors giving a new result stored in c - UNUSED in this program */ void add(float a[3], float b[3], float c[3]) { c[0] = a[0] + b[0] ; c[1] = a[1] + b[1] ; c[2] = a[2] + b[2] ; } /* subtract two vectors giving new vector c that points from b to a */ void subtract(float a[3], float b[3], float c[3]) { c[0] = a[0] - b[0] ; c[1] = a[1] - b[1] ; c[2] = a[2] - b[2] ; } /* dot product is equivalent to len(a)*len(b)*cos(anglebetweenab) - UNUSED*/ float dot(float a[3], float b[3]) { return( a[0]*b[0] + a[1]*b[1] + a[2]*b[2] ) ; } /* the cross product is a new vector c that is perpendicular to the a b plane */ void cross(float a[3], float b[3], float c[3]) { c[0] = a[1]*b[2] - a[2]*b[1] ; c[1] = a[2]*b[0] - a[0]*b[2] ; c[2] = a[0]*b[1] - a[1]*b[0] ; } void draw_string(float x, float y, float z, char *s) { glColor3f(0,1,1); glRasterPos3f(x, y, z) ; // place left baseline text start position while(*s) glutBitmapCharacter(GLUT_BITMAP_9_BY_15, *s++); } void draw_axes(void) { float k = 5 ; float nzero = 0.1 ; glDisable(GL_LIGHTING) ; glLineWidth(3) ; glBegin(GL_LINES) ; glColor3f(1,0,0) ; glVertex3f(nzero,nzero,nzero) ; glVertex3f(k,nzero,nzero) ; glColor3f(0,1,0) ; glVertex3f(nzero,nzero,nzero) ; glVertex3f(nzero,k,nzero) ; glColor3f(0,0,1) ; glVertex3f(nzero,nzero,nzero) ; glVertex3f(nzero,nzero,k) ; glEnd() ; glLineWidth(1) ; glColor3f(1,1,1) ; glRasterPos3f(k+1,0,0) ; glutBitmapCharacter(GLUT_BITMAP_HELVETICA_18,'X') ; glRasterPos3f(0,k+1,0) ; glutBitmapCharacter(GLUT_BITMAP_HELVETICA_18,'Y') ; glRasterPos3f(0,0,k+1) ; glutBitmapCharacter(GLUT_BITMAP_HELVETICA_18,'Z') ; glEnable(GL_LIGHTING) ; return ; } void draw_landscape(void) { int x,y ; float k ; glBegin(GL_QUADS) ; for( y = 0 ; y < HF_ROW ; y++ ) { for( x = 0 ; x < HF_COL ; x++ ) { glColor3fv(hf[y][x].color) ; glNormal3fv(hf[y][x].normal) ; glVertex3f(x,y, hf[y][x].height) ; glColor3fv(hf[y][x+1].color) ; glNormal3fv(hf[y][x+1].normal) ; glVertex3f(x+1,y, hf[y][x+1].height) ; glColor3fv(hf[y+1][x+1].color) ; glNormal3fv(hf[y+1][x+1].normal) ; glVertex3f(x+1,y+1, hf[y+1][x+1].height) ; glColor3fv(hf[y+1][x].color) ; glNormal3fv(hf[y+1][x].normal) ; glVertex3f(x,y+1, hf[y+1][x].height) ; } } glEnd() ; glEnable(GL_BLEND) ; glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA) ; k = sin((float)ticks/60.0) ; glBegin(GL_QUADS) ; glColor4f(0.0,0.0,0.5,0.5) ; glNormal3f(0,0,1) ; glVertex3f(0,0,k) ; glVertex3f(HF_COL,0,k) ; glVertex3f(HF_COL,HF_ROW,k) ; glVertex3f(0,HF_ROW,k) ; glEnd() ; glDisable(GL_BLEND) ; return ; } void setup_light(void) { float lightpos[4] = { 10, 10, 10, 0 } ; // w = 1 means the light glLightfv(GL_LIGHT0,GL_POSITION, lightpos) ; return ; } void draw(void) { float fromx, fromy, fromz ; float upx, upy, upz ; fromx = cam_r * sin(cam_phi) * cos(cam_theta) ; fromy = cam_r * sin(cam_phi) * sin(cam_theta) ; fromz = cam_r * cos(cam_phi) ; upx = sin(cam_phi+PI/2.0) * cos(cam_theta) ; upy = sin(cam_phi+PI/2.0) * sin(cam_theta) ; upz = cos(cam_phi+PI/2.0) ; glMatrixMode(GL_MODELVIEW); glLoadIdentity(); gluLookAt(fromx, fromy, fromz, 0,0,0, upx, upy, upz); setup_light() ; if( polymode == 0 ) glPolygonMode(GL_FRONT_AND_BACK, GL_FILL) ; else glPolygonMode(GL_FRONT_AND_BACK, GL_LINE) ; glPushMatrix() ; glTranslatef(-HF_COL/2, -HF_ROW/2, 0 ) ; draw_landscape() ; glPopMatrix() ; if( draw_axes_flag ) draw_axes() ; return ; } /* * GLUT Callback Stubs */ void display(void) { glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT); draw() ; glutSwapBuffers(); glutReportErrors(); } void timer(int value) { cam_theta += cam_dtheta ; if( (cam_phi+cam_dphi < 0) && (cam_phi+cam_dphi > -PI) ) cam_phi += cam_dphi ; else cam_dphi = 0 ; if( (cam_r+cam_dr < 800.0 ) && (cam_r+cam_dr > 1.0) ) cam_r += cam_dr ; else cam_dr = 0 ; if( go ) glutTimerFunc( 25, timer, value++); glutPostRedisplay() ; if( cam_dtheta ) cam_dtheta -= cam_dtheta/8.0 ; if( cam_dphi ) cam_dphi -= cam_dphi/8.0 ; if( cam_dr > 0 ) cam_dr -= 1.0/16.0 ; else if( cam_dr < 0 ) cam_dr += 1.0/16.0 ; ticks++ ; return ; } void reshape(int w, int h) { viewport_x = w ; viewport_y = h ; glViewport(0,0,viewport_x,viewport_y); glMatrixMode(GL_PROJECTION); glLoadIdentity(); gluPerspective(45,(float)w/(float)h,2,1000) ; return ; } void init_gl(int argc, char **argv) { glClearColor(0.25, 0.25, 0.25, 1); glEnable(GL_DEPTH_TEST); glEnable(GL_LIGHTING); glEnable(GL_LIGHT0); glEnable(GL_COLOR_MATERIAL); glShadeModel(GL_SMOOTH) ; } void apply_hts() { int x,y ; int side ; float k ; for(y = 0; y < HF_ROW; y++) { for(x = 0; x < HF_COL; x++) { float h = simg[y*ht_x+x]/6. - 2; hf[y][x].height = h; } } return; } char face_normal(int x, int y, float normal[3]) { float a[3],b[3],c[3] ; float v1[3],v2[3] ; if( (x < 0) || (y < 0) ) return(0) ; if( (x+1 > HF_COL) || (y+1 > HF_ROW) ) return(0) ; a[0] = x ; a[1] = y ; a[2] = hf[y][x].height ; b[0] = x+1 ; b[1] = y ; b[2] = hf[y][x+1].height ; c[0] = x+1 ; c[1] = y+1 ; c[2] = hf[y+1][x+1].height ; subtract(a,b,v1) ; subtract(a,c,v2) ; cross(v1,v2,normal) ; normalize(normal) ; return(1) ; } void generate_attribs(void) { int x,y ; float n[3],normal[3] ; float k ; for( y = 0 ; y < HF_ROW ; y++ ) { for( x = 0 ; x < HF_COL ; x++ ) { normal[0] = normal[1] = normal[2] = 0 ; k = 0 ; n[0] = n[1] = n[2] = 0 ; k += face_normal(x, y, n) ; normal[0] += n[0] ; normal[1] += n[1] ; normal[2] += n[2] ; n[0] = n[1] = n[2] = 0 ; k += face_normal(x-1, y, n) ; normal[0] += n[0] ; normal[1] += n[1] ; normal[2] += n[2] ; n[0] = n[1] = n[2] = 0 ; k += face_normal(x+1, y, n) ; normal[0] += n[0] ; normal[1] += n[1] ; normal[2] += n[2] ; n[0] = n[1] = n[2] = 0 ; k += face_normal(x, y-1, n) ; normal[0] += n[0] ; normal[1] += n[1] ; normal[2] += n[2] ; n[0] = n[1] = n[2] = 0 ; k += face_normal(x, y+1, n) ; normal[0] += n[0] ; normal[1] += n[1] ; normal[2] += n[2] ; normal[0] /= k ; normal[1] /= k ; normal[2] /= k ; hf[y][x].normal[0] = normal[0] ; hf[y][x].normal[1] = normal[1] ; hf[y][x].normal[2] = normal[2] ; k = hf[y][x].height + ((4.0*rand()/(RAND_MAX+1.0))-2.0) ; if( k > 50 ) { hf[y][x].color[0] = 1 ; hf[y][x].color[1] = 1 ; hf[y][x].color[2] = 1 ; } else if( k > 16 ) { float low[3],high[3] ; float t ; low[0] = 0 ; low[1] = 0.05 ; low[2] = 0 ; high[0] = 1 ; high[1] = 1 ; high[2] = 1 ; t = (k - 16.0) / (50.0-16.0) ; hf[y][x].color[0] = high[0]*t + (1.0-t)*low[0] ; hf[y][x].color[1] = high[1]*t + (1.0-t)*low[1] ; hf[y][x].color[2] = high[2]*t + (1.0-t)*low[2] ; } else if( k > 10 ) { float low[3],high[3] ; float t ; low[0] = 0 ; low[1] = 0.25 ; low[2] = 0 ; high[0] = 0 ; high[1] = 0.05 ; high[2] = 0 ; t = (k - 10.0) / (16.0-10.0) ; hf[y][x].color[0] = high[0]*t + (1.0-t)*low[0] ; hf[y][x].color[1] = high[1]*t + (1.0-t)*low[1] ; hf[y][x].color[2] = high[2]*t + (1.0-t)*low[2] ; } else if( k > 5 ) { float low[3],high[3] ; float t ; low[0] = 0.96 ; low[1] = 0.64 ; low[2] = 0.37 ; high[0] = 0 ; high[1] = 0.25 ; high[2] = 0 ; t = (k - 5.0) / (10.0-5.0) ; hf[y][x].color[0] = high[0]*t + (1.0-t)*low[0] ; hf[y][x].color[1] = high[1]*t + (1.0-t)*low[1] ; hf[y][x].color[2] = high[2]*t + (1.0-t)*low[2] ; } else if( k > 0 ) { float low[3],high[3] ; float t ; low[0] = 0.96 ; low[1] = 0.84 ; low[2] = 0.57 ; high[0] = 0.96 ; high[1] = 0.64 ; high[2] = 0.37 ; t = (k - 0.0) / (5.0-0.0) ; hf[y][x].color[0] = high[0]*t + (1.0-t)*low[0] ; hf[y][x].color[1] = high[1]*t + (1.0-t)*low[1] ; hf[y][x].color[2] = high[2]*t + (1.0-t)*low[2] ; } else if( k > -5 ) { float low[3],high[3] ; float t ; low[0] = 0 ; low[1] = 0 ; low[2] = 1 ; high[0] = 0.96 ; high[1] = 0.84 ; high[2] = 0.57 ; t = (k + 5.0) / (0.0+5.0) ; hf[y][x].color[0] = high[0]*t + (1.0-t)*low[0] ; hf[y][x].color[1] = high[1]*t + (1.0-t)*low[1] ; hf[y][x].color[2] = high[2]*t + (1.0-t)*low[2] ; } else { hf[y][x].color[0] = 0 ; hf[y][x].color[1] = 0 ; hf[y][x].color[2] = 1 ; } } } return ; } void refresh_view() { int i, n = ht_x*ht_y; for(i = 0; i < n; i++) // make an image copy as shorts for the wavelet simg[i] = ht_bytes[i]; update_hts(); // forward & reverse wavelet + your knockouts apply_hts(); // put the reconstructed hts on the landscape generate_attribs() ; // make the normals anad colors return ; } void keyboard( unsigned char key, int x, int y ) { switch( key ) { case 27: // escape key to exit exit(0); break; case ',': // rotate camera counter clockwise cam_dtheta -= cam_d/4.0 ; break ; case '.': // rotate camera clockwise cam_dtheta += cam_d/4.0 ; break ; case 'k': // tilt camera upwards cam_dphi += cam_d/4.0 ; break; case 'm': // tilt camera downwards cam_dphi -= cam_d/4.0 ; break; case '=': // move camera in toward middle cam_dr -= 0.25 ; break ; case '-': // move camera radially outward cam_dr += 0.25 ; break ; case 'a': // toggle axis drawing draw_axes_flag = ! draw_axes_flag ; break ; case 'p': // toggle pause go = !go ; fprintf(stderr,"%s\n", go?"run":"pause"); if( go ) glutTimerFunc( 25, timer, 0); break ; case 'l': // lower detail lev_detail /= 2; refresh_view(); break ; case 'h': // higher detail lev_detail += lev_detail;// double it if(lev_detail == 0) // special case recover from 0 lev_detail = 1; if(lev_detail > ht_x) // limit to full detail lev_detail = ht_x; refresh_view(); break ; case 'q': mode = ! mode ; break ; case 'w': polymode = ! polymode ; break ; default: //fprintf(stderr, "Unknown keystroke\n"); break; } glutPostRedisplay(); } void init(int argc, char **argv) { return ; } /* 1 = success, 0 = failure */ int pgm_p5_read(char *filename, int *x, int *y, unsigned char **bytes) { FILE *file = NULL ; char buf[1024] ; int colors, bpp ; int i ; file = fopen(filename,"rb") ; if( file == NULL ) { fprintf(stderr,"Unable to open '%s'\n",filename) ; return(0) ; } fscanf(file,"%s ",buf) ; if( strncmp("P5",buf,2) != 0 ) { fprintf(stderr,"P5 not found in '%s` header\n",filename) ; fclose(file) ; return(0) ; } /* ignoring the possible comments in the file here */ fscanf(file,"%d %d ",x, y) ; /* ignoring the possible comments in the file here */ fscanf(file,"%d",&colors) ; if( colors != 255 ) { fprintf(stderr,"In file '%s` header:\n",filename) ; fprintf(stderr,"Encountered a non-supported maximum color-component value '%d'\n",colors) ; fclose(file) ; return(0) ; } /* read the last '\n' */ fgetc(file) ; *bytes = (unsigned char *) malloc( (*x) * (*y) ) ; //for( i = (*y) - 1 ; i >= 0 ; i-- ) //fread(*bytes + i*(*x), (*x),1,file) ; fread(*bytes, (*x)*(*y), 1, file); fclose(file) ; return(1) ; } void init_postgl(void) { /* read in the image data */ pgm_p5_read("grandcan.pgm", &ht_x, &ht_y, &ht_bytes) ; /* allocate the array of shorts for the wavelet transform */ simg = (short *)malloc(ht_x*ht_y*sizeof(short)); lev_detail = ht_x; // start at perfect full detail reconstruction refresh_view(); } void init_glut(int argc, char **argv) { glutInit(&argc, argv); glutInitDisplayMode(GLUT_RGB|GLUT_DEPTH|GLUT_DOUBLE); glutInitWindowSize(WINDOW_X,WINDOW_Y); glutCreateWindow("Assignment 04"); glutDisplayFunc( display ); glutReshapeFunc( reshape ); glutKeyboardFunc( keyboard ); glutTimerFunc( 25, timer, 0); } int main(int argc, char **argv) { init(argc, argv); init_glut(argc, argv); init_gl(argc, argv); init_postgl() ; glutReportErrors(); glutMainLoop(); return(0) ; } #define OFF (262144+128) #define N 2048 // big enough for large ht map short lowdat[N/2+4]; // overshoot both ends short *low = lowdat+2; short high[N/2]; short pred[N/2]; fwstage(short *sig, int len, int step, int div) { int i; len /= 2; for(i = 0; i < len; i++) low[i] = (sig[(2*i)*step] + sig[(2*i+1)*step] + OFF)/2 - OFF/2; low[-1] = low[-2] = low[0]; low[len] = low[len+1] = low[len-1]; for(i = 0; i < len; i++) { pred[i] = 0; #define PRED // XXXX test setting this to NOPRED to degrade to Haar #ifdef PRED pred[i] += 14*low[i-2]; pred[i] -= 95*low[i-1]; pred[i] += 95*low[i+1]; pred[i] -= 14*low[i+2]; pred[i] /= 256; #endif } for(i = 0; i < len; i++) high[i] = (sig[(2*i+1)*step] - sig[(2*i)*step]) - pred[i]; for(i = 0; i < len; i++) { sig[i*step] = low[i]; sig[(i+len)*step] = high[i]; } } rwstage(short *sig, int len, int step, int div) { int i, v; len /= 2; for(i = 0; i < len; i++) { low[i] = sig[i*step]; high[i] = sig[(i+len)*step]; } low[-1] = low[-2] = low[0]; low[len] = low[len+1] = low[len-1]; for(i = 0; i < len; i++) { pred[i] = 0; #ifdef PRED pred[i] += 14*low[i-2]; pred[i] -= 95*low[i-1]; pred[i] += 95*low[i+1]; pred[i] -= 14*low[i+2]; pred[i] /= 256; #endif } for(i = 0; i < len; i++) { v = low[i] - ((high[i] + pred[i] + OFF)/2 - OFF/2); sig[(2*i)*step] = v; v = low[i] + ((high[i] + pred[i] + OFF + 1)/2 - OFF/2); if(div == 1 && step == 3) { if(v > 255) // fix over shoots v = 255; if(v < 0) v = 0; } sig[(2*i+1)*step] = v; } } update_hts() { int i, l, d, x, y; fprintf(stderr, "update_hts lev %d\n", lev_detail); // if you have transform problems just return from here /* forward transforms */ for(d = 1; d <= ht_x; d *= 2) { for(y = 0; y < ht_y/d; y++) fwstage(&simg[y*ht_x], ht_x/d, 1, d); for(x = 0; x < ht_x/d/2; x++) fwstage(&simg[x], ht_y/d, ht_x, d); for(x = ht_x/d/2; x < ht_x/d; x++) fwstage(&simg[x], ht_y/d, ht_x, d); } /* XXX do your knockout of high coefficients according to lev_detail */ /* if lev_detail is 0 then everything goes */ /* if it is, for example, 512 then keep only the top left 512*512 */ /* reverse transforms */ for(d = ht_x; d >= 1; d /= 2) { for(x = ht_x/d/2; x < ht_x/d; x++) rwstage(&simg[x], ht_y/d, ht_x, d); for(x = 0; x < ht_x/d/2; x++) rwstage(&simg[x], ht_y/d, ht_x, d); for(y = 0; y < ht_y/d; y++) rwstage(&simg[y*ht_x], ht_x/d, 1, d); } }