Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

glutil.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: David Woolford, 11/06/2007 (woolford@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #ifdef EMAN2_USING_OPENGL
00037 
00038 #ifdef _WIN32
00039         #include <windows.h>
00040 #endif
00041 
00042 #ifndef GL_GLEXT_PROTOTYPES
00043         #define GL_GLEXT_PROTOTYPES
00044 #endif  //GL_GLEXT_PROTOTYPES
00045 
00046 #include "glutil.h"
00047 #include "emdata.h"
00048 #include "marchingcubes.h"
00049 
00050 #ifdef __APPLE__
00051         #include "OpenGL/gl.h"
00052         #include "OpenGL/glu.h"
00053         #include "OpenGL/glext.h"
00054 #else // WIN32, LINUX
00055         #include "GL/gl.h"
00056         #include "GL/glu.h"
00057         #include "GL/glext.h"
00058 #endif  //__APPLE__
00059 
00060 using namespace EMAN;
00061 
00062 //By depfault we need to first bind data to the GPU
00063 GLuint GLUtil::buffer[2] = {0, 0};
00064 
00065 unsigned int GLUtil::gen_glu_mipmaps(const EMData* const emdata)
00066 {
00067         if ( emdata->get_data() == 0 ) throw NullPointerException("Error, attempt to create an OpenGL mipmap without internally stored data");
00068         ENTERFUNC;
00069 
00070         unsigned int tex_name;
00071         glGenTextures(1, &tex_name);
00072 
00073         if ( emdata->ny == 1 && emdata->nz == 1 )
00074         {
00075                 glBindTexture(GL_TEXTURE_1D, tex_name);
00076                 gluBuild1DMipmaps(GL_TEXTURE_1D, GL_LUMINANCE, emdata->nx, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00077         } else if (emdata->nz == 1) {
00078                 glBindTexture(GL_TEXTURE_2D, tex_name);
00079                 gluBuild2DMipmaps(GL_TEXTURE_2D, GL_LUMINANCE, emdata->nx, emdata->ny, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00080         }
00081         else {
00082 #ifdef  _WIN32
00083                 //There is no gluBuild3DMipmaps() function in glu.h for VS2003 and VS2005
00084                 printf("3D OpenGL mipmaps are not available on this platform.\n");
00085 #else
00086                 glBindTexture(GL_TEXTURE_3D, tex_name);
00087                 gluBuild3DMipmaps(GL_TEXTURE_3D, GL_LUMINANCE, emdata->nx, emdata->ny, emdata->nz, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00088 #endif  //_WIN32
00089         }
00090 
00091         EXITFUNC;
00092         return tex_name;
00093 }
00094 
00095 unsigned int GLUtil::gen_gl_texture(const EMData* const emdata, GLenum format)
00096 {
00097         if ( emdata->get_data() == 0 ) throw NullPointerException("Error, attempt to create an OpenGL texture without internally stored data");
00098         ENTERFUNC;
00099 
00100         unsigned int tex_name;
00101         glGenTextures(1, &tex_name);
00102 
00103         if ( emdata->ny == 1 && emdata->nz == 1 )
00104         {
00105                 glBindTexture(GL_TEXTURE_1D, tex_name);
00106                 glTexImage1D(GL_TEXTURE_1D,0,format,emdata->nx,0,format, GL_FLOAT, (void*)(emdata->get_data()));
00107         } else if (emdata->nz == 1) {
00108                 glBindTexture(GL_TEXTURE_2D, tex_name);
00109                 glTexImage2D(GL_TEXTURE_2D,0,format,emdata->nx,emdata->ny,0,format, GL_FLOAT, (void*)(emdata->get_data()));
00110         }
00111         else {
00112                 glBindTexture(GL_TEXTURE_3D, tex_name);
00113 #ifdef _WIN32
00114         PFNGLTEXIMAGE3DPROC glTexImage3D;
00115 #endif
00116                 glTexImage3D(GL_TEXTURE_3D,0, format,emdata->nx,emdata->ny,emdata->nz,0,format, GL_FLOAT, (void*)(emdata->get_data()));
00117         }
00118 
00119         EXITFUNC;
00120         return tex_name;
00121 }
00122 
00123 unsigned int GLUtil::render_amp8_gl_texture(EMData* emdata, int x0, int y0, int ixsize, int iysize, int bpl, float scale, int mingray, int maxgray,     float render_min, float render_max,float gamma,int flags) {
00124 
00125         string pixels = render_amp8(emdata, x0, y0, ixsize,iysize, bpl, scale, mingray, maxgray, render_min, render_max, gamma,flags);
00126 
00127         unsigned int tex_name;
00128         glGenTextures(1, &tex_name);
00129 
00130         glBindTexture(GL_TEXTURE_2D,tex_name);
00131         glTexImage2D(GL_TEXTURE_2D,0,GL_LUMINANCE,ixsize,iysize,0,GL_LUMINANCE, GL_UNSIGNED_BYTE, pixels.c_str() );
00132 
00133         return tex_name;
00134 }
00135 
00136 // undef GL_GLEXT_PROTOTYPES
00137 #ifdef GL_GLEXT_PROTOTYPES
00138 #undef GL_GLEXT_PROTOTYPES
00139 #endif
00140 
00141 int GLUtil::nearest_projected_points(const vector<float>& model_matrix, const vector<float>& proj_matrix, const vector<int>& view_matrix,
00142                 const vector<Vec3f>& points, const float mouse_x, const float mouse_y,const float& nearness)
00143 {
00144         double proj[16];
00145         double model[16];
00146         int view[4];
00147 
00148         copy(proj_matrix.begin(), proj_matrix.end(), proj);
00149         copy(model_matrix.begin(), model_matrix.end(), model);
00150         copy(view_matrix.begin(), view_matrix.end(), view);
00151 
00152         vector<Vec3f> unproj_points;
00153         double x,y,z;
00154         double r,s,t;
00155         for(vector<Vec3f>::const_iterator it = points.begin(); it != points.end(); ++it) {
00156                 r = (double) (*it)[0];
00157                 s = (double) (*it)[1];
00158                 t = (double) (*it)[2];
00159                 gluProject(r,s,t, model, proj, view, &x,&y,&z);
00160                 unproj_points.push_back(Vec3f(x,y,z));
00161         }
00162 
00163         vector<int> intersections;
00164 
00165         float n_squared = nearness*nearness;
00166         for(unsigned int i = 0; i < unproj_points.size(); ++i){
00167                 Vec3f& v = unproj_points[i];
00168                 float dx = v[0] - mouse_x;
00169                 dx *= dx;
00170                 float dy = v[1] - mouse_y;
00171                 dy *= dy;
00172                 if ((dx+dy) <= n_squared) intersections.push_back((int)i);
00173         }
00174 
00175         int intersection = -1;
00176         float near_z = 0;
00177         for(vector<int>::const_iterator it = intersections.begin(); it != intersections.end(); ++it) {
00178                 if (intersection == -1 || unproj_points[*it][2] < near_z) {
00179                         intersection = *it;
00180                         near_z = unproj_points[*it][2];
00181                 }
00182         }
00183 
00184         return intersection;
00185 }
00186 
00187 void GLUtil::colored_rectangle(const vector<float>& data,const float& alpha,const bool center_point){
00188 
00189         glBegin(GL_LINE_LOOP);
00190         glColor4f(data[0],data[1],data[2],alpha);
00191         glVertex2i(data[3],data[4]);
00192         glVertex2i(data[5],data[4]);
00193         glVertex2i(data[5],data[6]);
00194         glVertex2i(data[3],data[6]);
00195         glEnd();
00196 
00197         if (center_point) {
00198                 glBegin(GL_POINTS);
00199                 glVertex2f( (data[3]+data[5])/2, (data[4]+data[6])/2);
00200                 glEnd();
00201         }
00202 }
00203 
00204 void GLUtil::mx_bbox(const vector<float>& data, const vector<float>& text_color, const vector<float>& bg_color) {
00205         glDisable(GL_TEXTURE_2D);
00206         glBegin(GL_QUADS);
00207         if (bg_color.size()==4) glColor4f(bg_color[0],bg_color[1],bg_color[2],bg_color[3]);
00208         else glColor3f(bg_color[0],bg_color[1],bg_color[2]);
00209         glVertex3f(data[0]-1,data[1]-1,-.1f);
00210         glVertex3f(data[3]+1,data[1]-1,-.1f);
00211         glVertex3f(data[3]+1,data[4]+1,-.1f);
00212         glVertex3f(data[0]-1,data[4]+1,-.1f);
00213         glEnd();
00214 //      glEnable(GL_TEXTURE_2D);
00215         if (text_color.size()==4) glColor4f(text_color[0],text_color[1],text_color[2],text_color[3]);
00216         else glColor3f(text_color[0],text_color[1],text_color[2]);
00217         
00218 }
00219 
00220 std::string GLUtil::render_amp8(EMData* emdata, int x0, int y0, int ixsize, int iysize,
00221                                                  int bpl, float scale, int mingray, int maxgray,
00222                                                  float render_min, float render_max,float gamma,int flags)
00223 {
00224         ENTERFUNC;
00225 
00226         int asrgb;
00227         int hist=(flags&2)/2;
00228         int invy=(flags&4)?1:0;
00229 
00230         int nx = emdata->nx;
00231         int ny = emdata->ny;
00232 //      int nz = emdata->nz;
00233         int nxy = emdata->nx * emdata->ny;
00234 
00235         if (emdata->get_ndim() > 2) {
00236                 throw ImageDimensionException("1D/2D only");
00237         }
00238 
00239         if (emdata->is_complex()) {
00240                 emdata->ri2ap();
00241         }
00242 
00243         if (render_max <= render_min) {
00244                 render_max = render_min + 0.01f;
00245         }
00246 
00247         if (gamma<=0) gamma=1.0;
00248 
00249         // Calculating a full floating point gamma for
00250         // each piGLUtil::xel in the image slows rendering unacceptably
00251         // however, applying a gamma-mapping to an 8 bit colorspace
00252         // has unacceptable coarse accuracy. So, we oversample the 8 bit colorspace
00253         // as a 12 bit colorspace and apply the gamma mapping to that
00254         // This should produce good accuracy for gamma values
00255         // larger than 0.5 (and a high upper limit)
00256         static int smg0=0,smg1=0;       // while this destroys threadsafety in the rendering process
00257         static float sgam=0;            // it is necessary for speed when rendering large numbers of small images
00258         static unsigned char gammamap[4096];
00259         if (gamma!=1.0 && (smg0!=mingray || smg1!=maxgray || sgam!=gamma)) {
00260                 for (int i=0; i<4096; i++) {
00261                         if (mingray<maxgray) gammamap[i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00262                         else gammamap[4095-i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00263                 }
00264         }
00265         smg0=mingray;   // so we don't recompute the map unless something changes
00266         smg1=maxgray;
00267         sgam=gamma;
00268 
00269         if (flags&8) asrgb=4;
00270         else if (flags&1) asrgb=3;
00271         else asrgb=1;
00272 
00273         std::string ret=std::string();
00274 //      ret.resize(iysize*bpl);
00275         ret.assign(iysize*bpl+hist*1024,char(mingray));
00276         unsigned char *data=(unsigned char *)ret.data();
00277         unsigned int *histd=(unsigned int *)(data+iysize*bpl);
00278         if (hist) {
00279                 for (int i=0; i<256; i++) histd[i]=0;
00280         }
00281 
00282         float rm = render_min;
00283         float inv_scale = 1.0f / scale;
00284         int ysize = iysize;
00285         int xsize = ixsize;
00286 
00287         int ymin = 0;
00288         if (iysize * inv_scale > ny) {
00289                 ymin = (int) (iysize - ny / inv_scale);
00290         }
00291 
00292         float gs = (maxgray - mingray) / (render_max - render_min);
00293         float gs2 = 4095.999f / (render_max - render_min);
00294 //      float gs2 = 1.0 / (render_max - render_min);
00295         if (render_max < render_min) {
00296                 gs = 0;
00297                 rm = FLT_MAX;
00298         }
00299 
00300         int dsx = -1;
00301         int dsy = 0;
00302         int remx = 0;
00303         int remy = 0;
00304         const int scale_n = 100000;
00305 
00306         int addi = 0;
00307         int addr = 0;
00308         if (inv_scale == floor(inv_scale)) {
00309                 dsx = (int) inv_scale;
00310                 dsy = (int) (inv_scale * nx);
00311         }
00312         else {
00313                 addi = (int) floor(inv_scale);
00314                 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00315         }
00316 
00317         int xmin = 0;
00318         if (x0 < 0) {
00319                 xmin = (int) (-x0 / inv_scale);
00320                 xsize -= (int) floor(x0 / inv_scale);
00321                 x0 = 0;
00322         }
00323 
00324         if ((xsize - xmin) * inv_scale > (nx - x0)) {
00325                 xsize = (int) ((nx - x0) / inv_scale + xmin);
00326         }
00327         int ymax = ysize - 1;
00328         if (y0 < 0) {
00329                 ymax = (int) (ysize + y0 / inv_scale - 1);
00330                 ymin += (int) floor(y0 / inv_scale);
00331                 y0 = 0;
00332         }
00333 
00334         if (xmin < 0) xmin = 0;
00335         if (ymin < 0) ymin = 0;
00336         if (xsize > ixsize) xsize = ixsize;
00337         if (ymax > iysize) ymax = iysize;
00338 
00339         int lmax = nx * ny - 1;
00340 
00341         int mid=nx*ny/2;
00342         float * image_data = emdata->get_data();
00343 
00345 
00346         const int rangemax=4082;
00347         int gaussianwide=5000;
00348         unsigned int* grayhe=NULL;
00349         float* gaussianpdf=NULL;
00350         float* gaussiancdf=NULL;
00351         int* gaussianlookup=NULL;
00352         
00353         unsigned int* graypdftwo=NULL;
00354         
00355         bool binflag=0;
00356         int binlocation=-1;
00357         
00358         if (flags&32){
00359                 int graypdf[rangemax]={0};//256
00360                 int graycdf[rangemax-2]={0};//254
00361                 //unsigned int grayhe[rangemax-2]={0};//#number=254
00362                 
00363                 graypdftwo = new unsigned int[maxgray-mingray];
00364                 grayhe = new unsigned int[rangemax-2];
00365                 //unsigned char graylookup[(int)(render_max-render_min)];//render_max-render_min
00366                 
00367                 for (int i=0; i<(int)(rangemax-2); i++) {
00368                         grayhe[i] = 0;    // Initialize all elements to zero.
00369                 }
00370                 
00371                 for (int i=0; i<(maxgray-mingray); i++) {//0~253
00372                         graypdftwo[i]=0;
00373                 }
00374                 
00375                 if (dsx != -1) {
00376                         int l = x0 + y0 * nx;
00377                         for (int j = ymax; j >= ymin; j--) {
00378                                 int br = l;
00379                                 for (int i = xmin; i < xsize; i++) {
00380                                         if (l > lmax) {
00381                                                 break;
00382                                         }
00383                                         int k = 0;
00384                                         int p;
00385                                         float t;
00386                                         if (dsx==1) t=image_data[l];
00387                                         else {                                          // This block does local pixel averaging for nicer reduced views
00388                                                 t=0;
00389                                                 for (int iii=0; iii<dsx; iii++) {
00390                                                         for (int jjj=0; jjj<dsy; jjj+=nx) {
00391                                                                 t+=image_data[l+iii+jjj];
00392                                                         }
00393                                                 }
00394                                                 t/=dsx*(dsy/nx);
00395                                         }
00396 
00397                                         if (t <= rm) graypdf[0]++;
00398                                         else if (t >= render_max) graypdf[rangemax-1]++;                        
00399                                         else {
00400                                                 graypdf[(int)(ceil((rangemax-2)*(t - render_min)/(render_max-render_min)))]++;
00401                                                 graypdftwo[(unsigned char) (gs * (t - render_min))]++;
00402                                         }
00403                                         l += dsx;
00404                                 }
00405                                 l = br + dsy;
00406                         }
00407                 }
00408                 else {
00409                         remy = 10;
00410                         int l = x0 + y0 * nx;
00411                         for (int j = ymax; j >= ymin; j--) {
00412                                 int addj = addi;
00413                                 // There seems to be some overflow issue happening
00414                                 // where the statement if (l > lmax) break (below) doesn't work
00415                                 // because the loop that iterates jjj can inadvertantly go out of bounds
00416                                 if (( l + addi*nx ) >= nxy ) {
00417                                         addj = (nxy-l)/nx;
00418                                         if (addj <= 0) continue;
00419                                 }
00420                                 int br = l;
00421                                 remx = 10;
00422                                 for (int i = xmin; i < xsize; i++) {
00423                                         if (l > lmax) break;
00424                                         int k = 0;
00425                                         int p;
00426                                         float t;
00427                                         if (addi<=1) t = image_data[l];
00428                                         else {                                          // This block does local pixel averaging for nicer reduced views
00429                                                 t=0;
00430                                                 for (int jjj=0; jjj<addj; jjj++) {
00431                                                         for (int iii=0; iii<addi; iii++) {
00432                                                                 t+=image_data[l+iii+jjj*nx];
00433                                                         }
00434                                                 }
00435                                                 t/=addi*addi;
00436                                         }
00438                                         if (t <= rm) graypdf[0]++;
00439                                         else if (t >= render_max) graypdf[rangemax-1]++;
00440                                         else {
00441                                                 graypdf[(int)(ceil((rangemax-2)*(t - render_min)/(render_max-render_min)))]++;
00442                                         }
00443 
00444                                         data[i * asrgb + j * bpl] = p;
00445                                         if (hist) histd[p]++;
00446                                         l += addi;
00447                                         remx += addr;
00448                                         if (remx > scale_n) {
00449                                                 remx -= scale_n;
00450                                                 l++;
00451                                         }
00452                                 }
00453                                 l = br + addi * nx;
00454                                 remy += addr;
00455                                 if (remy > scale_n) {
00456                                         remy -= scale_n;
00457                                         l += nx; 
00458                                 }
00459                         }
00460                 }
00461                 
00462                 for (int i=0; i<(rangemax-2); i++) {//0~253
00463                         for (int j=0;j<(i+1);j++) {
00464                                 graycdf[i]=graycdf[i]+graypdf[j+1];
00465                         }
00466                 }
00467                 
00468                 //graypdftwo binflag
00469                 binflag=0;
00470                 for (int i=1; i<(maxgray-mingray); i++) {//0~253
00471                         
00472                         if (((float)graypdftwo[i]/graycdf[rangemax-3])>0.2){
00473                                 binflag=1;
00474                                 binlocation=i;
00475                                 break;
00476                         }
00477                 }               
00478                 
00479                 if (binflag==1){
00480                         for (int i=(binlocation*16+1); i<((binlocation+1)*16); i++) {
00481                                 graypdf[i]=0;
00482                                 //graypdf[i]=(graycdf[rangemax-3]-graypdftwo[binlocation])/(maxgray-mingray);
00483                         }
00484 
00485                         for (int i=0; i<(rangemax-2); i++) {//0~253
00486                                 graycdf[i]=0;
00487                         }
00488                         
00489                         for (int i=0; i<(rangemax-2); i++) {//0~253
00490                                 for (int j=0;j<(i+1);j++) {
00491                                         graycdf[i]=graycdf[i]+graypdf[j+1];
00492                                 }
00493                         }
00494                 }
00495                 
00496                 
00497                 // start gaussian matching
00498                 float mean=abs(rangemax-2)/2;
00499                 float standdv=abs(mean)/3;
00500                 gaussianpdf=new float[rangemax-2];
00501                 gaussiancdf=new float[rangemax-2];
00502                 gaussianlookup=new int[rangemax-2];
00503                 
00504                 for (int i=0; i<(rangemax-2); i++){
00505                         gaussianpdf[i]=exp(-(i-mean)*(i-mean)/(2*standdv*standdv))/sqrt(standdv * standdv * 2 * M_PI);
00506                         
00507                         if(i!=0){
00508                                 gaussiancdf[i]=gaussiancdf[i-1]+gaussianpdf[i];
00509                         }
00510                         else{
00511                                 gaussiancdf[i]=gaussianpdf[i];
00512                         }
00513                 }
00514                 
00515                 for (int i=0; i<(rangemax-2); i++) {
00516                         gaussiancdf[i]=graycdf[rangemax-3]*gaussiancdf[i]/gaussiancdf[rangemax-3];
00517                 }
00518                 
00519                 for (int i=0; i<(rangemax-2); i++){
00520                         for (int j=0; j<(rangemax-2); j++){
00521                                 if (graycdf[i]<=gaussiancdf[j]){
00522                                         gaussianlookup[i]=j;
00523                                         break;
00524                                 }
00525                         }
00526                 }
00527                 
00528                 for (int i=0; i<(rangemax-2); i++) {                    
00529                         grayhe[i]=floor(0.5+(((double)(rangemax-3)*graycdf[i])/graycdf[rangemax-3]));
00530                 }
00531                 
00532         } 
00533                 
00535         
00536         
00537         if (emdata->is_complex()) {
00538                 if (dsx != -1) {
00539                         int l = y0 * nx;
00540                         for (int j = ymax; j >= ymin; j--) {
00541                                 int ll = x0;
00542                                 for (int i = xmin; i < xsize; i++) {
00543                                         if (l + ll > lmax || ll >= nx - 2) break;
00544 
00545                                         int k = 0;
00546                                         unsigned char p;
00547                                         if (ll >= nx / 2) {
00548                                                 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
00549                                                 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00550                                         }
00551                                         else k = nx * ny - (l + 2 * ll) - 2;
00552                                         if (k>=mid) k-=mid;             // These 2 lines handle the Fourier origin being in the corner, not the middle
00553                                         else k+=mid;
00554                                         float t = image_data[k];
00555                                         int ph;
00556                                         // in color mode
00557                                         if (flags&16 && asrgb>2) {
00558 //                                              if (l >= (ny - inv_scale) * nx) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;  // complex phase as integer 0-767;
00559                                                 if (ll >= nx / 2) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;        // complex phase as integer 0-767;
00560                                                 else ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;   // complex phase as integer 0-767;
00561                                         }
00562                                         if (t <= rm)  p = mingray;
00563                                         else if (t >= render_max) p = maxgray;
00564                                         else if (gamma!=1.0) {
00565                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00566                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00567 //                                              p = (unsigned char) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00568 //                                              p += mingray;
00569 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00570 //                                              k += mingray;
00571                                         }
00572                                         else {
00573                                                 p = (unsigned char) (gs * (t - render_min));
00574                                                 p += mingray;
00575                                         }
00576                                         // color rendering
00577                                         if (flags&16 && asrgb>2) {
00578                                                 if (ph<256) {
00579                                                         data[i * asrgb + j * bpl] = p*(255-ph)/256;
00580                                                         data[i * asrgb + j * bpl+1] = p*ph/256;
00581                                                         data[i * asrgb + j * bpl+2] = 0;
00582                                                 }
00583                                                 else if (ph<512) {
00584                                                         data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00585                                                         data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00586                                                         data[i * asrgb + j * bpl] = 0;
00587                                                 }
00588                                                 else {
00589                                                         data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00590                                                         data[i * asrgb + j * bpl] = p*(ph-512)/256;
00591                                                         data[i * asrgb + j * bpl+1] = 0;
00592                                                 }
00593                                         }
00594                                         else data[i * asrgb + j * bpl] = p;
00595                                         if (hist) histd[p]++;
00596                                         ll += dsx;
00597                                 }
00598                                 l += dsy;
00599                         }
00600                 }
00601                 else {
00602                         remy = 10;
00603                         int l = y0 * nx;
00604                         for (int j = ymax; j >= ymin; j--) {
00605                                 int br = l;
00606                                 remx = 10;
00607                                 int ll = x0;
00608                                 for (int i = xmin; i < xsize - 1; i++) {
00609                                         if (l + ll > lmax || ll >= nx - 2) {
00610                                                 break;
00611                                         }
00612                                         int k = 0;
00613                                         unsigned char p;
00614                                         if (ll >= nx / 2) {
00615                                                 if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
00616                                                 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00617                                         }
00618                                         else k = nx * ny - (l + 2 * ll) - 2;
00619                                         if (k>=mid) k-=mid;             // These 2 lines handle the Fourier origin being in the corner, not the middle
00620                                         else k+=mid;
00621 
00622                                         float t = image_data[k];
00623                                         // in color mode
00624                                         int ph;
00625                                         if (flags&16 && asrgb>2) {
00626                                                 if (l >= (ny * nx - nx)) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;         // complex phase as integer 0-767;
00627                                                 else ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;   // complex phase as integer 0-767;
00628                                         }
00629                                         if (t <= rm)
00630                                                 p = mingray;
00631                                         else if (t >= render_max) {
00632                                                 p = maxgray;
00633                                         }
00634                                         else if (gamma!=1.0) {
00635                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00636                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00637 //                                              p = (unsigned char) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00638 //                                              p += mingray;
00639 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00640 //                                              k += mingray;
00641                                         }
00642                                         else {
00643                                                 p = (unsigned char) (gs * (t - render_min));
00644                                                 p += mingray;
00645                                         }
00646                                         if (flags&16 && asrgb>2) {
00647                                                 if (ph<256) {
00648                                                         data[i * asrgb + j * bpl] = p*(255-ph)/256;
00649                                                         data[i * asrgb + j * bpl+1] = p*ph/256;
00650                                                         data[i * asrgb + j * bpl+2] = 0;
00651                                                 }
00652                                                 else if (ph<512) {
00653                                                         data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00654                                                         data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00655                                                         data[i * asrgb + j * bpl] = 0;
00656                                                 }
00657                                                 else {
00658                                                         data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00659                                                         data[i * asrgb + j * bpl] = p*(ph-512)/256;
00660                                                         data[i * asrgb + j * bpl+1] = 0;
00661                                                 }
00662                                         }
00663                                         else data[i * asrgb + j * bpl] = p;
00664                                         if (hist) histd[p]++;
00665                                         ll += addi;
00666                                         remx += addr;
00667                                         if (remx > scale_n) {
00668                                                 remx -= scale_n;
00669                                                 ll++;
00670                                         }
00671                                 }
00672                                 l = br + addi * nx;
00673                                 remy += addr;
00674                                 if (remy > scale_n) {
00675                                         remy -= scale_n;
00676                                         l += nx;
00677                                 }
00678                         }
00679                 }
00680         }
00681         else {
00682                 if (dsx != -1) {
00683                         int l = x0 + y0 * nx;
00684                         for (int j = ymax; j >= ymin; j--) {
00685                                 int br = l;
00686                                 for (int i = xmin; i < xsize; i++) {
00687                                         if (l > lmax) {
00688                                                 break;
00689                                         }
00690                                         int k = 0;
00691                                         unsigned char p;
00692                                         float t;
00693                                         if (dsx==1) t=image_data[l];
00694                                         else {                                          // This block does local pixel averaging for nicer reduced views
00695                                                 t=0;
00696                                                 for (int iii=0; iii<dsx; iii++) {
00697                                                         for (int jjj=0; jjj<dsy; jjj+=nx) {
00698                                                                 t+=image_data[l+iii+jjj];
00699                                                         }
00700                                                 }
00701                                                 t/=dsx*(dsy/nx);
00702                                         }
00703 
00704                                         if (t <= rm) p = mingray;
00705                                         else if (t >= render_max) p = maxgray;
00706                                         else if (gamma!=1.0) {
00707                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00708                                                 p = gammamap[k];// apply gamma using precomputed gamma map
00709 //                                              k = (int) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00710 //                                              k += mingray;
00711 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00712 //                                              k += mingray;
00713                                         }
00714                                         else {
00715                                                 if (flags&32){
00716                                                         //p = graylookup[(int)(t - render_min)];
00717                                                         //graylookup[i]=(unsigned char)((maxgray-mingray-2)*grayhe[(int)(i*(rangemax-3)/(render_max-render_min))]/(rangemax-3)+1);
00718                                                         if (flags&64)
00719                                                         {
00720                                                                 p=(unsigned char)(gaussianlookup[(int)(floor((rangemax-2)*(t - render_min)/(render_max-render_min)))]*(maxgray-mingray-2)/(rangemax-3)+1);
00721                                                         }
00722                                                         else{
00723                                                                 p=(unsigned char)(grayhe[(int)((t - render_min)*(rangemax-3)/(render_max-render_min))]*(maxgray-mingray-2)/(rangemax-3)+1);}
00724                                                         //p=(unsigned char)gaussianlookup[(int)(ceil)((t - render_min)*(rangemax-3)/(render_max-render_min))]+1;
00725                                                 }
00726                                                 else
00727                                                 {
00728                                                         p=(unsigned char) (gs * (t - render_min));
00729                                                 }
00730                                                 
00731                                                 p += mingray;
00732                                         }
00733                                         data[i * asrgb + j * bpl] = p;
00734                                         if (hist) histd[p]++;
00735                                         l += dsx;
00736                                 }
00737                                 l = br + dsy;
00738                         }
00739                 }
00740                 else {
00741                         remy = 10;
00742                         int l = x0 + y0 * nx;
00743                         for (int j = ymax; j >= ymin; j--) {
00744                                 int addj = addi;
00745                                 // There seems to be some overflow issue happening
00746                                 // where the statement if (l > lmax) break (below) doesn't work
00747                                 // because the loop that iterates jjj can inadvertantly go out of bounds
00748                                 if (( l + addi*nx ) >= nxy ) {
00749                                         addj = (nxy-l)/nx;
00750                                         if (addj <= 0) continue;
00751                                 }
00752                                 int br = l;
00753                                 remx = 10;
00754                                 for (int i = xmin; i < xsize; i++) {
00755                                         if (l > lmax) break;
00756                                         int k = 0;
00757                                         unsigned char p;
00758                                         float t;
00759                                         if (addi<=1) t = image_data[l];
00760                                         else {                                          // This block does local pixel averaging for nicer reduced views
00761                                                 t=0;
00762                                                 for (int jjj=0; jjj<addj; jjj++) {
00763                                                         for (int iii=0; iii<addi; iii++) {
00764                                                                 t+=image_data[l+iii+jjj*nx];
00765                                                         }
00766                                                 }
00767                                                 t/=addi*addi;
00768                                         }
00769                                         if (t <= rm) p = mingray;
00770                                         else if (t >= render_max) p = maxgray;
00771                                         else if (gamma!=1.0) {
00772                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00773                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00774 //                                              k = (int) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00775 //                                              k += mingray;
00776 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00777 //                                              k += mingray;
00778                                         }
00779                                         else {
00780                                                 if (flags&32){
00781                                                         //p = graylookup[(int)(t - render_min)];
00782                                                         if (flags&64){
00783                                                                 p=(unsigned char)(gaussianlookup[(int)(floor((rangemax-2)*(t - render_min)/(render_max-render_min)))]*(maxgray-mingray-2)/(rangemax-3)+1);
00784                                                         }
00785                                                         else{
00786                                                                 p=(unsigned char)(grayhe[(int)((t - render_min)*(rangemax-3)/(render_max-render_min))]*(maxgray-mingray-2)/(rangemax-3)+1);
00787                                                         }
00788                                                         //p=(unsigned char)gaussianlookup[(int)(ceil)((t - render_min)*(rangemax-3)/(render_max-render_min))]+1;
00789                                                 }
00790                                                 else{
00791                                                         p=(unsigned char) (gs * (t - render_min));      
00792                                                 }
00793                                                 
00794                                                 p += mingray;
00795                                         }
00796                                         data[i * asrgb + j * bpl] = p;
00797                                         if (hist) histd[p]++;
00798                                         l += addi;
00799                                         remx += addr;
00800                                         if (remx > scale_n) {
00801                                                 remx -= scale_n;
00802                                                 l++;
00803                                         }
00804                                 }
00805                                 l = br + addi * nx;
00806                                 remy += addr;
00807                                 if (remy > scale_n) {
00808                                         remy -= scale_n;
00809                                         l += nx;
00810                                 }
00811                         }
00812                 }
00813         }
00814 
00815         // this replicates r -> g,b
00816         if (asrgb==3 && !(flags&16)) {
00817                 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00818                         for (int i=xmin; i<xsize*3; i+=3) {
00819                                 data[i+j+1]=data[i+j+2]=data[i+j];
00820                         }
00821                 }
00822         }
00823         if (asrgb==4 && !(flags&16)) {
00824                 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00825                         for (int i=xmin; i<xsize*4; i+=4) {
00826                                 data[i+j+1]=data[i+j+2]=data[i+j+3]=data[i+j];
00827                                 data[i+j+3]=255;
00828                         }
00829                 }
00830         }
00831 
00832         EXITFUNC;
00833 
00834         // ok, ok, not the most efficient place to do this, but it works
00835         if (invy) {
00836                 int x,y;
00837                 char swp;
00838                 for (y=0; y<iysize/2; y++) {
00839                         for (x=0; x<ixsize; x++) {
00840                                 swp=ret[y*bpl+x];
00841                                 ret[y*bpl+x]=ret[(iysize-y-1)*bpl+x];
00842                                 ret[(iysize-y-1)*bpl+x]=swp;
00843                         }
00844                 }
00845         }
00846 
00847     //  return PyString_FromStringAndSize((const char*) data,iysize*bpl);
00848         if (flags&16) {
00849                 glDrawPixels(ixsize,iysize,GL_LUMINANCE,GL_UNSIGNED_BYTE,(const GLvoid *)ret.data());
00850         }
00851 
00852 
00853         delete[] grayhe;
00854         delete[] gaussianpdf;
00855         delete[] gaussiancdf;
00856         delete[] gaussianlookup;
00857         delete[] graypdftwo;
00858         
00859         return ret;
00860 }
00861 
00862 //DEPRECATED
00863 unsigned long GLUtil::get_isosurface_dl(MarchingCubes* mc, unsigned int tex_id,bool surface_face_z, bool recontour)
00864 {
00865         if ( mc->_isodl != 0 ) glDeleteLists(mc->_isodl,1);
00866                 
00867         if ( recontour ) mc->calculate_surface();
00868         if ( surface_face_z ) mc->surface_face_z();
00869         if( mc->getRGBmode() ) mc->color_vertices();
00870         
00871 #if MARCHING_CUBES_DEBUG
00872         cout << "There are " << ff.elem()/3 << " faces and " << pp.elem() << " points and " << nn.elem() << " normals to render in generate dl" << endl;
00873 #endif
00874         int maxf;
00875 //#ifdef        _WIN32
00876 //      glGetIntegerv(GL_MAX_ELEMENTS_INDICES_WIN,&maxf);
00877 //#else
00878         glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
00879 //#endif        //_WIN32
00880 #if MARCHING_CUBES_DEBUG
00881         int maxv;
00882         glGetIntegerv(GL_MAX_ELEMENTS_VERTICES,&maxv);
00883         cout << "Max vertices is " << maxv << " max indices is " << maxf << endl;
00884         cout << "Using OpenGL " << glGetString(GL_VERSION) << endl;
00885 #endif
00886 
00887         if ( recontour ) for (unsigned int i = 0; i < mc->ff.elem(); ++i ) mc->ff[i] /= 3;
00888 
00889         if ( maxf % 3 != 0 )
00890         {
00891                 maxf = maxf - (maxf%3);
00892         }
00893 
00894         if ( tex_id != 0 ) {
00895                 // Normalize the coordinates to be on the interval 0,1
00896                 mc->pp.mult3(1.0f/(float) mc->_emdata->get_xsize(), 1.0f/(float)mc->_emdata->get_ysize(), 1.0f/(float)mc->_emdata->get_zsize());
00897                 glEnableClientState(GL_TEXTURE_COORD_ARRAY);
00898                 glDisableClientState(GL_NORMAL_ARRAY);
00899                 glTexCoordPointer(3,GL_FLOAT,0,mc->pp.get_data());
00900         }
00901         else {
00902                 glEnableClientState(GL_NORMAL_ARRAY);
00903                 glDisableClientState(GL_TEXTURE_COORD_ARRAY);
00904                 glNormalPointer(GL_FLOAT,0, mc->nn.get_data());
00905         }
00906 
00907         glEnableClientState(GL_VERTEX_ARRAY);
00908         glVertexPointer(3,GL_FLOAT,0, mc->pp.get_data());
00909         
00910         if (mc->getRGBmode()){
00911                 glEnableClientState(GL_COLOR_ARRAY);
00912                 glColorPointer(3,GL_FLOAT,0, mc->cc.get_data());
00913         }
00914         
00915         mc->_isodl = glGenLists(1);
00916         
00917 #if MARCHING_CUBES_DEBUG
00918         int time0 = clock();
00919 #endif
00920         
00921         glNewList(mc->_isodl,GL_COMPILE);
00922 
00923         if ( tex_id != 0 ) {
00924                 glEnable(GL_TEXTURE_3D);
00925                 glBindTexture(GL_TEXTURE_3D, tex_id);
00926                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);
00927                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);
00928                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP);
00929                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
00930                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
00931                 glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
00932         }
00933         
00934         // Drawing range elements based on the output of glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
00935         // Saved about 60% of the time... drawRange should probably always be true
00936         bool drawRange = true;
00937         if ( drawRange == false ) {
00938                 glDrawElements(GL_TRIANGLES,mc->ff.elem(),GL_UNSIGNED_INT, mc->ff.get_data());
00939         } else {
00940                 for(unsigned int i = 0; i < mc->ff.elem(); i+=maxf)
00941                 {
00942                         if ( (i+maxf) > mc->ff.elem())
00943                                 glDrawElements(GL_TRIANGLES,mc->ff.elem()-i,GL_UNSIGNED_INT,&(mc->ff[i]));
00944                         else
00945                                 glDrawElements(GL_TRIANGLES,maxf,GL_UNSIGNED_INT,&(mc->ff[i]));
00946                                 // glDrawRangeElements is part of the extensions, we might want to experiment with its performance at some stage,
00947                                 // so please leave this code here, commented out. This is an either or situation, so if glDrawRangeElements is used,
00948                                 // glDrawElements above would have to be commented out.
00949                                 // glDrawRangeElements(GL_TRIANGLES,0,0,maxf,GL_UNSIGNED_INT,&ff[i]);
00950                 }
00951         }
00952 
00953         if ( tex_id != 0 ) glDisable(GL_TEXTURE_3D);
00954 
00955         glEndList();
00956         
00957 #if MARCHING_CUBES_DEBUG
00958         int time1 = clock();
00959         cout << "It took " << (time1-time0) << " " << (float)(time1-time0)/CLOCKS_PER_SEC << " to draw elements" << endl;
00960 #endif
00961         return mc->_isodl;
00962 }
00963 
00964 void GLUtil::contour_isosurface(MarchingCubes* mc)
00965 {
00966         mc->calculate_surface();
00967         //What does this do???
00968         for (unsigned int i = 0; i < mc->ff.elem(); ++i ) mc->ff[i] /= 3;
00969         //Need to rebind data (to the GPU)
00970         mc->needtobind = true;
00971 }
00972 
00973 void GLUtil::render_using_VBOs(MarchingCubes* mc, unsigned int tex_id,bool surface_face_z)
00974 {
00975         // In current version Texture is not supported b/c it is not used... EVER
00976         
00977 #ifdef _WIN32
00978         typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
00979         PFNGLGENBUFFERSPROC glGenBuffers;
00980         glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
00981 
00982         typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
00983         PFNGLISBUFFERPROC glIsBuffer;
00984         glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
00985 #endif  //_WIN32
00986 
00987         if ( surface_face_z ) mc->surface_face_z();
00988         
00989         //Bug in OpenGL, sometimes glGenBuffers doesn't work. Try again, if we still fail then return...
00990         if (!glIsBuffer(mc->buffer[0])) glGenBuffers(4, mc->buffer);
00991 
00992         //whenever something changes, like color mode or color scale (or threshold), we need to recolor
00993         if( mc->getRGBmode() && (mc->rgbgenerator.getNeedToRecolor() || mc->needtobind)){
00994                 mc->color_vertices();
00995                 mc->needtobind = true;
00996         }
00997 
00998         int maxf;
00999         glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
01000 
01001         if ( maxf % 3 != 0 )
01002         {
01003                 maxf = maxf - (maxf%3);
01004         }
01005 
01006 #ifdef _WIN32
01007         typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
01008         PFNGLBINDBUFFERPROC glBindBuffer;
01009         glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
01010 
01011         typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size, const GLvoid *data, GLenum usage);
01012         PFNGLBUFFERDATAPROC glBufferData;
01013         glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
01014 #endif  //_WIN32
01015 
01016         //Normal vectors
01017         glEnableClientState(GL_NORMAL_ARRAY);
01018         glDisableClientState(GL_TEXTURE_COORD_ARRAY);
01019         glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[2]);
01020         if (mc->needtobind) glBufferData(GL_ARRAY_BUFFER, 4*mc->nn.elem(), mc->nn.get_data(), GL_STATIC_DRAW);
01021         glNormalPointer(GL_FLOAT,0,0);
01022 
01023         //Vertex vectors
01024         glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[0]);
01025         if (mc->needtobind) glBufferData(GL_ARRAY_BUFFER, 4*mc->pp.elem(), mc->pp.get_data(), GL_STATIC_DRAW);
01026         glEnableClientState(GL_VERTEX_ARRAY);
01027         glVertexPointer(3,GL_FLOAT,0,0);
01028 
01029         //Color vectors
01030         if (mc->getRGBmode()){
01031                 glEnableClientState(GL_COLOR_ARRAY);
01032                 glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[3]);
01033                 if (mc->needtobind) glBufferData(GL_ARRAY_BUFFER, 4*mc->cc.elem(), mc->cc.get_data(), GL_STATIC_DRAW);
01034                 glColorPointer(3,GL_FLOAT,0, 0);
01035         }
01036         
01037         //Indices
01038         glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, mc->buffer[1]);
01039         if (mc->needtobind) glBufferData(GL_ELEMENT_ARRAY_BUFFER, 4*mc->ff.elem(), mc->ff.get_data(), GL_STATIC_DRAW);
01040         
01041         // This lets us know if buffers an not implmemted
01042         if (!glIsBuffer(mc->buffer[0])){
01043                 cout << "Can't Generate GL Vertex Buffers. glGenBuffer error" << endl;
01044                 return;
01045         }
01046         
01047         //finally draw the elemenets
01048         glDrawElements(GL_TRIANGLES,mc->ff.elem(),GL_UNSIGNED_INT,0);
01049         // No longer need to bind data to the GPU
01050         mc->needtobind = false;
01051 
01052 }
01053 
01054 // This crap could be avoided and speed up(A lot) if we implemented an openGL shader........
01055 void GLUtil::glLoadMatrix(const Transform& xform)
01056 {
01057 #ifdef _WIN32
01058         typedef void (APIENTRYP PFNGLLOADTRANSPOSEMATRIXFPROC) (const GLfloat *m);
01059         PFNGLLOADTRANSPOSEMATRIXFPROC glLoadTransposeMatrixf = NULL;
01060         glLoadTransposeMatrixf = (PFNGLLOADTRANSPOSEMATRIXFPROC) wglGetProcAddress("glLoadTransposeMatrixf");
01061 #endif  //_WIN32
01062 
01063         vector<float> xformlist = xform.get_matrix_4x4();
01064         glLoadTransposeMatrixf(reinterpret_cast<GLfloat*>(&xformlist[0]));
01065 }
01066 
01067 void GLUtil::glMultMatrix(const Transform& xform)
01068 {
01069 #ifdef _WIN32
01070         typedef void (APIENTRYP PFNGLMULTTRANSPOSEMATRIXFPROC) (const GLfloat *m);
01071         PFNGLMULTTRANSPOSEMATRIXFPROC glMultTransposeMatrixf = NULL;
01072         glMultTransposeMatrixf = (PFNGLMULTTRANSPOSEMATRIXFPROC) wglGetProcAddress("glMultTransposeMatrixf");
01073 #endif  //_WIN32
01074 
01075         vector<float> xformlist = xform.get_matrix_4x4();
01076         glMultTransposeMatrixf(reinterpret_cast<GLfloat*>(&xformlist[0]));
01077 }
01078 
01079 //This draw a bounding box, or any box 
01080 void GLUtil::glDrawBoundingBox(float width, float height, float depth)
01081 {
01082         
01083         float w2 = width/2.0f;
01084         float h2 = height/2.0f;
01085         float d2 = depth/2.0f;
01086         
01087         float vertices[24] = {-w2,  h2,  d2, w2,  h2,  d2, w2, -h2,  d2, -w2, -h2,  d2, -w2,  h2, -d2, w2,  h2, -d2, w2, -h2, -d2, -w2, -h2, -d2};
01088         int indices[24] = {0, 1, 1, 2, 2, 3, 3, 0, 4, 5, 5, 6, 6, 7, 7, 4, 0, 4, 3, 7, 1, 5, 2, 6};
01089         
01090 #ifdef _WIN32
01091         typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
01092         PFNGLGENBUFFERSPROC glGenBuffers;
01093         glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
01094 
01095         typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
01096         PFNGLISBUFFERPROC glIsBuffer;
01097         glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
01098 
01099         typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
01100         PFNGLBINDBUFFERPROC glBindBuffer;
01101         glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
01102 
01103         typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size, const GLvoid *data, GLenum usage);
01104         PFNGLBUFFERDATAPROC glBufferData;
01105         glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
01106 #endif  //_WIN32
01107 
01108         if ( glIsBuffer(GLUtil::buffer[0])  == GL_FALSE ){
01109                 glGenBuffers(2, GLUtil::buffer);
01110         }
01111         
01112         // Could use dirty bit here but not worth my time to implment
01113         glBindBuffer(GL_ARRAY_BUFFER, GLUtil::buffer[0]);
01114         glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);
01115         glEnableClientState(GL_VERTEX_ARRAY);
01116         glVertexPointer(3,GL_FLOAT,0,0);
01117         
01118         glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, GLUtil::buffer[1]);
01119         glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(indices), indices, GL_STATIC_DRAW);
01120         glDrawElements(GL_LINES,24,GL_UNSIGNED_INT,0);
01121         
01122 }
01123 
01124 void GLUtil::glDrawDisk(float radius, int spokes)
01125 {
01126 #ifdef _WIN32
01127         typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
01128         PFNGLGENBUFFERSPROC glGenBuffers;
01129         glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
01130 
01131         typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
01132         PFNGLISBUFFERPROC glIsBuffer;
01133         glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
01134 
01135         typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
01136         PFNGLBINDBUFFERPROC glBindBuffer;
01137         glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
01138 
01139         typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size, const GLvoid *data, GLenum usage);
01140         PFNGLBUFFERDATAPROC glBufferData;
01141         glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
01142 #endif  //_WIN32
01143 
01144         //This code is experimental
01145         int arraysize = 4*pow((double)2,(double)spokes);
01146         int sideofarray = pow((double)2,(double)spokes);
01147         
01148 //      float vertices[3*arraysize + 3];
01149         vector<float> vertices(3*arraysize + 3);
01150         
01151         //last vertex is center
01152         vertices[3*arraysize] = 0.0;
01153         vertices[3*arraysize+1] = 0.0;
01154         vertices[3*arraysize+2] = 0.0;
01155         
01156         vertices[0] = 0.0;
01157         vertices[1] = radius;
01158         vertices[2] = 0.0;
01159                 
01160         vertices[sideofarray*3] = radius;
01161         vertices[sideofarray*3 + 1] = 0.0;
01162         vertices[sideofarray*3 + 2] = 0.0;
01163         
01164         vertices[sideofarray*6] = 0.0;
01165         vertices[sideofarray*6 + 1] = -radius;
01166         vertices[sideofarray*6 + 2] = 0.0;
01167         
01168         vertices[sideofarray*9] = -radius;
01169         vertices[sideofarray*9 + 1] = 0.0;
01170         vertices[sideofarray*9 + 2] = 0.0;
01171         
01172         //This could aslo be implemented recusivly
01173         for(int step = 0; step < spokes; step++){
01174                 //starting location
01175                 int x = sideofarray/pow(2.0,(double)(step+1));
01176                 for(int i = 1; i <= 4*pow(2.0,(double)step); i++ ){
01177                         
01178                         // take the necessary steps
01179                         int index =  x + 2*x*(i-1);
01180                         int index_f = (index + x) % arraysize;
01181                         int index_i = index - x;
01182                         cout << index << " " << index_f << " " << index_i << endl;
01183                         
01184                         //need to resclae length to that of radius
01185                         vertices[index_f*3] = (vertices[index_f*3] - vertices[index_i*3])/2.0f;
01186                         vertices[index_f*3 + 1] = (vertices[index_f*3 + 1] - vertices[index_i*3 + 1])/2.0f;
01187                         vertices[index_f*3 + 2] = (vertices[index_f*3 + 2] - vertices[index_i*3 + 2])/2.0f;
01188                 }
01189         }
01190         
01191         //GL stuff
01192         if ( glIsBuffer(GLUtil::buffer[0])  == GL_FALSE ){
01193                 glGenBuffers(2, GLUtil::buffer);
01194         }
01195         
01196         // Could use dirty bit here but not worth my time to implment
01197         glBindBuffer(GL_ARRAY_BUFFER, GLUtil::buffer[0]);
01198         //glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);
01199         glBufferData(GL_ARRAY_BUFFER, vertices.size(), &(vertices[0]), GL_STATIC_DRAW);
01200         glEnableClientState(GL_VERTEX_ARRAY);
01201         glVertexPointer(3,GL_FLOAT,0,0);
01202         
01203         //Code to select indices
01204         
01205         
01206 }
01207 #endif // EMAN2_USING_OPENGL

Generated on Tue Jun 11 13:40:38 2013 for EMAN2 by  doxygen 1.3.9.1