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 #include "glutil.h"
00039 #include "emdata.h"
00040 #include "marchingcubes.h"
00041 
00042 #ifndef _WIN32
00043         #ifndef GL_GLEXT_PROTOTYPES
00044                 #define GL_GLEXT_PROTOTYPES
00045         #endif  //GL_GLEXT_PROTOTYPES
00046 #endif  //_WIN32
00047 
00048 #ifdef __APPLE__
00049         #include "OpenGL/gl.h"
00050         #include "OpenGL/glu.h"
00051         #include "OpenGL/glext.h"
00052 #else // WIN32, LINUX
00053         #include "GL/gl.h"
00054         #include "GL/glu.h"
00055         #include "GL/glext.h"
00056 #endif  //__APPLE__
00057 
00058 using namespace EMAN;
00059 
00060 unsigned int GLUtil::gen_glu_mipmaps(const EMData* const emdata)
00061 {
00062         if ( emdata->get_data() == 0 ) throw NullPointerException("Error, attempt to create an OpenGL mipmap without internally stored data");
00063         ENTERFUNC;
00064 
00065         unsigned int tex_name;
00066         glGenTextures(1, &tex_name);
00067 
00068         if ( emdata->ny == 1 && emdata->nz == 1 )
00069         {
00070                 glBindTexture(GL_TEXTURE_1D, tex_name);
00071                 gluBuild1DMipmaps(GL_TEXTURE_1D, GL_LUMINANCE, emdata->nx, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00072         } else if (emdata->nz == 1) {
00073                 glBindTexture(GL_TEXTURE_2D, tex_name);
00074                 gluBuild2DMipmaps(GL_TEXTURE_2D, GL_LUMINANCE, emdata->nx, emdata->ny, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00075         }
00076         else {
00077 #ifdef  _WIN32
00078                 //There is no gluBuild3DMipmaps() function in glu.h for VS2003 and VS2005
00079                 printf("3D OpenGL mipmaps are not available on this platform.\n");
00080 #else
00081                 glBindTexture(GL_TEXTURE_3D, tex_name);
00082                 gluBuild3DMipmaps(GL_TEXTURE_3D, GL_LUMINANCE, emdata->nx, emdata->ny, emdata->nz, GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00083 #endif  //_WIN32
00084         }
00085 
00086         EXITFUNC;
00087         return tex_name;
00088 }
00089 
00090 unsigned int GLUtil::gen_gl_texture(const EMData* const emdata)
00091 {
00092         if ( emdata->get_data() == 0 ) throw NullPointerException("Error, attempt to create an OpenGL texture without internally stored data");
00093         ENTERFUNC;
00094 
00095         unsigned int tex_name;
00096         glGenTextures(1, &tex_name);
00097 
00098         if ( emdata->ny == 1 && emdata->nz == 1 )
00099         {
00100                 glBindTexture(GL_TEXTURE_1D, tex_name);
00101                 glTexImage1D(GL_TEXTURE_1D,0,GL_LUMINANCE,emdata->nx,0,GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00102         } else if (emdata->nz == 1) {
00103                 glBindTexture(GL_TEXTURE_2D, tex_name);
00104                 glTexImage2D(GL_TEXTURE_2D,0,GL_LUMINANCE,emdata->nx,emdata->ny,0,GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00105         }
00106         else {
00107                 glBindTexture(GL_TEXTURE_3D, tex_name);
00108 #ifdef _WIN32
00109         PFNGLTEXIMAGE3DPROC glTexImage3D;
00110 #endif
00111                 glTexImage3D(GL_TEXTURE_3D,0, GL_LUMINANCE,emdata->nx,emdata->ny,emdata->nz,0,GL_LUMINANCE, GL_FLOAT, (void*)(emdata->get_data()));
00112         }
00113 
00114         EXITFUNC;
00115         return tex_name;
00116 }
00117 
00118 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) {
00119 
00120         string pixels = render_amp8(emdata, x0, y0, ixsize,iysize, bpl, scale, mingray, maxgray, render_min, render_max, gamma,flags);
00121 
00122         unsigned int tex_name;
00123         glGenTextures(1, &tex_name);
00124 
00125         glBindTexture(GL_TEXTURE_2D,tex_name);
00126         glTexImage2D(GL_TEXTURE_2D,0,GL_LUMINANCE,ixsize,iysize,0,GL_LUMINANCE, GL_UNSIGNED_BYTE, pixels.c_str() );
00127 
00128         return tex_name;
00129 }
00130 
00131 // undef GL_GLEXT_PROTOTYPES
00132 #ifdef GL_GLEXT_PROTOTYPES
00133 #undef GL_GLEXT_PROTOTYPES
00134 #endif
00135 
00136 int GLUtil::nearest_projected_points(const vector<float>& model_matrix, const vector<float>& proj_matrix, const vector<int>& view_matrix,
00137                 const vector<Vec3f>& points, const float mouse_x, const float mouse_y,const float& nearness)
00138 {
00139         double proj[16];
00140         double model[16];
00141         int view[4];
00142 
00143         copy(proj_matrix.begin(), proj_matrix.end(), proj);
00144         copy(model_matrix.begin(), model_matrix.end(), model);
00145         copy(view_matrix.begin(), view_matrix.end(), view);
00146 
00147         vector<Vec3f> unproj_points;
00148         double x,y,z;
00149         double r,s,t;
00150         for(vector<Vec3f>::const_iterator it = points.begin(); it != points.end(); ++it) {
00151                 r = (double) (*it)[0];
00152                 s = (double) (*it)[1];
00153                 t = (double) (*it)[2];
00154                 gluProject(r,s,t, model, proj, view, &x,&y,&z);
00155                 unproj_points.push_back(Vec3f(x,y,z));
00156         }
00157 
00158         vector<int> intersections;
00159 
00160         float n_squared = nearness*nearness;
00161         for(unsigned int i = 0; i < unproj_points.size(); ++i){
00162                 Vec3f& v = unproj_points[i];
00163                 float dx = v[0] - mouse_x;
00164                 dx *= dx;
00165                 float dy = v[1] - mouse_y;
00166                 dy *= dy;
00167                 if ((dx+dy) <= n_squared) intersections.push_back((int)i);
00168         }
00169 
00170         int intersection = -1;
00171         float near_z = 0;
00172         for(vector<int>::const_iterator it = intersections.begin(); it != intersections.end(); ++it) {
00173                 if (intersection == -1 || unproj_points[*it][2] < near_z) {
00174                         intersection = *it;
00175                         near_z = unproj_points[*it][2];
00176                 }
00177         }
00178 
00179         return intersection;
00180 }
00181 
00182 void GLUtil::colored_rectangle(const vector<float>& data,const float& alpha,const bool center_point){
00183 
00184         glBegin(GL_LINE_LOOP);
00185         glColor4f(data[0],data[1],data[2],alpha);
00186         glVertex2i(data[3],data[4]);
00187         glVertex2i(data[5],data[4]);
00188         glVertex2i(data[5],data[6]);
00189         glVertex2i(data[3],data[6]);
00190         glEnd();
00191 
00192         if (center_point) {
00193                 glBegin(GL_POINTS);
00194                 glVertex2f( (data[3]+data[5])/2, (data[4]+data[6])/2);
00195                 glEnd();
00196         }
00197 }
00198 
00199 void GLUtil::mx_bbox(const vector<float>& data, const vector<float>& text_color, const vector<float>& bg_color) {
00200         glColor4f(bg_color[0],bg_color[1],bg_color[2],bg_color[3]);
00201         glDisable(GL_TEXTURE_2D);
00202         glBegin(GL_QUADS);
00203 
00204         glVertex3f(data[0]-1,data[1]-1,-.1);
00205         glVertex3f(data[3]+1,data[1]-1,-.1);
00206         glVertex3f(data[3]+1,data[4]+1,-.1);
00207         glVertex3f(data[0]-1,data[4]+1,-.1);
00208         glEnd();
00209         glEnable(GL_TEXTURE_2D);
00210         glColor4f(text_color[0],text_color[1],text_color[2],text_color[3]);
00211 }
00212 
00213 std::string GLUtil::render_amp8(EMData* emdata, int x0, int y0, int ixsize, int iysize,
00214                                                  int bpl, float scale, int mingray, int maxgray,
00215                                                  float render_min, float render_max,float gamma,int flags)
00216 {
00217         ENTERFUNC;
00218 
00219         int asrgb;
00220         int hist=(flags&2)/2;
00221         int invy=(flags&4)?1:0;
00222 
00223         int nx = emdata->nx;
00224         int ny = emdata->ny;
00225 //      int nz = emdata->nz;
00226         int nxy = emdata->nx * emdata->ny;
00227 
00228         if (emdata->get_ndim() > 2) {
00229                 throw ImageDimensionException("1D/2D only");
00230         }
00231 
00232         if (emdata->is_complex()) {
00233                 emdata->ri2ap();
00234         }
00235 
00236         if (render_max <= render_min) {
00237                 render_max = render_min + 0.01f;
00238         }
00239 
00240         if (gamma<=0) gamma=1.0;
00241 
00242         // Calculating a full floating point gamma for
00243         // each piGLUtil::xel in the image slows rendering unacceptably
00244         // however, applying a gamma-mapping to an 8 bit colorspace
00245         // has unacceptable coarse accuracy. So, we oversample the 8 bit colorspace
00246         // as a 12 bit colorspace and apply the gamma mapping to that
00247         // This should produce good accuracy for gamma values
00248         // larger than 0.5 (and a high upper limit)
00249         static int smg0=0,smg1=0;       // while this destroys threadsafety in the rendering process
00250         static float sgam=0;            // it is necessary for speed when rendering large numbers of small images
00251         static unsigned char gammamap[4096];
00252         if (gamma!=1.0 && (smg0!=mingray || smg1!=maxgray || sgam!=gamma)) {
00253                 for (int i=0; i<4096; i++) {
00254                         if (mingray<maxgray) gammamap[i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00255                         else gammamap[4095-i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00256                 }
00257         }
00258         smg0=mingray;   // so we don't recompute the map unless something changes
00259         smg1=maxgray;
00260         sgam=gamma;
00261 
00262         if (flags&8) asrgb=4;
00263         else if (flags&1) asrgb=3;
00264         else asrgb=1;
00265 
00266         std::string ret=std::string();
00267 //      ret.resize(iysize*bpl);
00268         ret.assign(iysize*bpl+hist*1024,char(mingray));
00269         unsigned char *data=(unsigned char *)ret.data();
00270         unsigned int *histd=(unsigned int *)(data+iysize*bpl);
00271         if (hist) {
00272                 for (int i=0; i<256; i++) histd[i]=0;
00273         }
00274 
00275         float rm = render_min;
00276         float inv_scale = 1.0f / scale;
00277         int ysize = iysize;
00278         int xsize = ixsize;
00279 
00280         int ymin = 0;
00281         if (iysize * inv_scale > ny) {
00282                 ymin = (int) (iysize - ny / inv_scale);
00283         }
00284 
00285         float gs = (maxgray - mingray) / (render_max - render_min);
00286         float gs2 = 4095.999f / (render_max - render_min);
00287 //      float gs2 = 1.0 / (render_max - render_min);
00288         if (render_max < render_min) {
00289                 gs = 0;
00290                 rm = FLT_MAX;
00291         }
00292 
00293         int dsx = -1;
00294         int dsy = 0;
00295         int remx = 0;
00296         int remy = 0;
00297         const int scale_n = 100000;
00298 
00299         int addi = 0;
00300         int addr = 0;
00301         if (inv_scale == floor(inv_scale)) {
00302                 dsx = (int) inv_scale;
00303                 dsy = (int) (inv_scale * nx);
00304         }
00305         else {
00306                 addi = (int) floor(inv_scale);
00307                 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00308         }
00309 
00310         int xmin = 0;
00311         if (x0 < 0) {
00312                 xmin = (int) (-x0 / inv_scale);
00313                 xsize -= (int) floor(x0 / inv_scale);
00314                 x0 = 0;
00315         }
00316 
00317         if ((xsize - xmin) * inv_scale > (nx - x0)) {
00318                 xsize = (int) ((nx - x0) / inv_scale + xmin);
00319         }
00320         int ymax = ysize - 1;
00321         if (y0 < 0) {
00322                 ymax = (int) (ysize + y0 / inv_scale - 1);
00323                 ymin += (int) floor(y0 / inv_scale);
00324                 y0 = 0;
00325         }
00326 
00327         if (xmin < 0) xmin = 0;
00328         if (ymin < 0) ymin = 0;
00329         if (xsize > ixsize) xsize = ixsize;
00330         if (ymax > iysize) ymax = iysize;
00331 
00332         int lmax = nx * ny - 1;
00333 
00334         int mid=nx*ny/2;
00335         float * image_data = emdata->get_data();
00336         if (emdata->is_complex()) {
00337                 if (dsx != -1) {
00338                         int l = y0 * nx;
00339                         for (int j = ymax; j >= ymin; j--) {
00340                                 int ll = x0;
00341                                 for (int i = xmin; i < xsize; i++) {
00342                                         if (l + ll > lmax || ll >= nx - 2) break;
00343 
00344                                         int k = 0;
00345                                         unsigned char p;
00346                                         if (ll >= nx / 2) {
00347                                                 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
00348                                                 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00349                                         }
00350                                         else k = nx * ny - (l + 2 * ll) - 2;
00351                                         if (k>=mid) k-=mid;             // These 2 lines handle the Fourier origin being in the corner, not the middle
00352                                         else k+=mid;
00353                                         float t = image_data[k];
00354                                         int ph;
00355                                         // in color mode
00356                                         if (flags&16 && asrgb>2) {
00357 //                                              if (l >= (ny - inv_scale) * nx) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;  // complex phase as integer 0-767;
00358                                                 if (ll >= nx / 2) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;        // complex phase as integer 0-767;
00359                                                 else ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;   // complex phase as integer 0-767;
00360                                         }
00361                                         if (t <= rm)  p = mingray;
00362                                         else if (t >= render_max) p = maxgray;
00363                                         else if (gamma!=1.0) {
00364                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00365                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00366 //                                              p = (unsigned char) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00367 //                                              p += mingray;
00368 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00369 //                                              k += mingray;
00370                                         }
00371                                         else {
00372                                                 p = (unsigned char) (gs * (t - render_min));
00373                                                 p += mingray;
00374                                         }
00375                                         // color rendering
00376                                         if (flags&16 && asrgb>2) {
00377                                                 if (ph<256) {
00378                                                         data[i * asrgb + j * bpl] = p*(255-ph)/256;
00379                                                         data[i * asrgb + j * bpl+1] = p*ph/256;
00380                                                         data[i * asrgb + j * bpl+2] = 0;
00381                                                 }
00382                                                 else if (ph<512) {
00383                                                         data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00384                                                         data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00385                                                         data[i * asrgb + j * bpl] = 0;
00386                                                 }
00387                                                 else {
00388                                                         data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00389                                                         data[i * asrgb + j * bpl] = p*(ph-512)/256;
00390                                                         data[i * asrgb + j * bpl+1] = 0;
00391                                                 }
00392                                         }
00393                                         else data[i * asrgb + j * bpl] = p;
00394                                         if (hist) histd[p]++;
00395                                         ll += dsx;
00396                                 }
00397                                 l += dsy;
00398                         }
00399                 }
00400                 else {
00401                         remy = 10;
00402                         int l = y0 * nx;
00403                         for (int j = ymax; j >= ymin; j--) {
00404                                 int br = l;
00405                                 remx = 10;
00406                                 int ll = x0;
00407                                 for (int i = xmin; i < xsize - 1; i++) {
00408                                         if (l + ll > lmax || ll >= nx - 2) {
00409                                                 break;
00410                                         }
00411                                         int k = 0;
00412                                         unsigned char p;
00413                                         if (ll >= nx / 2) {
00414                                                 if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
00415                                                 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00416                                         }
00417                                         else k = nx * ny - (l + 2 * ll) - 2;
00418                                         if (k>=mid) k-=mid;             // These 2 lines handle the Fourier origin being in the corner, not the middle
00419                                         else k+=mid;
00420 
00421                                         float t = image_data[k];
00422                                         // in color mode
00423                                         int ph;
00424                                         if (flags&16 && asrgb>2) {
00425                                                 if (l >= (ny * nx - nx)) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;         // complex phase as integer 0-767;
00426                                                 else ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;   // complex phase as integer 0-767;
00427                                         }
00428                                         if (t <= rm)
00429                                                 p = mingray;
00430                                         else if (t >= render_max) {
00431                                                 p = maxgray;
00432                                         }
00433                                         else if (gamma!=1.0) {
00434                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00435                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00436 //                                              p = (unsigned char) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00437 //                                              p += mingray;
00438 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00439 //                                              k += mingray;
00440                                         }
00441                                         else {
00442                                                 p = (unsigned char) (gs * (t - render_min));
00443                                                 p += mingray;
00444                                         }
00445                                         if (flags&16 && asrgb>2) {
00446                                                 if (ph<256) {
00447                                                         data[i * asrgb + j * bpl] = p*(255-ph)/256;
00448                                                         data[i * asrgb + j * bpl+1] = p*ph/256;
00449                                                         data[i * asrgb + j * bpl+2] = 0;
00450                                                 }
00451                                                 else if (ph<512) {
00452                                                         data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00453                                                         data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00454                                                         data[i * asrgb + j * bpl] = 0;
00455                                                 }
00456                                                 else {
00457                                                         data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00458                                                         data[i * asrgb + j * bpl] = p*(ph-512)/256;
00459                                                         data[i * asrgb + j * bpl+1] = 0;
00460                                                 }
00461                                         }
00462                                         else data[i * asrgb + j * bpl] = p;
00463                                         if (hist) histd[p]++;
00464                                         ll += addi;
00465                                         remx += addr;
00466                                         if (remx > scale_n) {
00467                                                 remx -= scale_n;
00468                                                 ll++;
00469                                         }
00470                                 }
00471                                 l = br + addi * nx;
00472                                 remy += addr;
00473                                 if (remy > scale_n) {
00474                                         remy -= scale_n;
00475                                         l += nx;
00476                                 }
00477                         }
00478                 }
00479         }
00480         else {
00481                 if (dsx != -1) {
00482                         int l = x0 + y0 * nx;
00483                         for (int j = ymax; j >= ymin; j--) {
00484                                 int br = l;
00485                                 for (int i = xmin; i < xsize; i++) {
00486                                         if (l > lmax) {
00487                                                 break;
00488                                         }
00489                                         int k = 0;
00490                                         unsigned char p;
00491                                         float t;
00492                                         if (dsx==1) t=image_data[l];
00493                                         else {                                          // This block does local pixel averaging for nicer reduced views
00494                                                 t=0;
00495                                                 for (int iii=0; iii<dsx; iii++) {
00496                                                         for (int jjj=0; jjj<dsy; jjj+=nx) {
00497                                                                 t+=image_data[l+iii+jjj];
00498                                                         }
00499                                                 }
00500                                                 t/=dsx*(dsy/nx);
00501                                         }
00502 
00503                                         if (t <= rm) p = mingray;
00504                                         else if (t >= render_max) p = maxgray;
00505                                         else if (gamma!=1.0) {
00506                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00507                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00508 //                                              k = (int) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00509 //                                              k += mingray;
00510 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00511 //                                              k += mingray;
00512                                         }
00513                                         else {
00514                                                 p = (unsigned char) (gs * (t - render_min));
00515                                                 p += mingray;
00516                                         }
00517                                         data[i * asrgb + j * bpl] = p;
00518                                         if (hist) histd[p]++;
00519                                         l += dsx;
00520                                 }
00521                                 l = br + dsy;
00522                         }
00523                 }
00524                 else {
00525                         remy = 10;
00526                         int l = x0 + y0 * nx;
00527                         for (int j = ymax; j >= ymin; j--) {
00528                                 int addj = addi;
00529                                 // There seems to be some overflow issue happening
00530                                 // where the statement if (l > lmax) break (below) doesn't work
00531                                 // because the loop that iterates jjj can inadvertantly go out of bounds
00532                                 if (( l + addi*nx ) >= nxy ) {
00533                                         addj = (nxy-l)/nx;
00534                                         if (addj <= 0) continue;
00535                                 }
00536                                 int br = l;
00537                                 remx = 10;
00538                                 for (int i = xmin; i < xsize; i++) {
00539                                         if (l > lmax) break;
00540                                         int k = 0;
00541                                         unsigned char p;
00542                                         float t;
00543                                         if (addi<=1) t = image_data[l];
00544                                         else {                                          // This block does local pixel averaging for nicer reduced views
00545                                                 t=0;
00546                                                 for (int jjj=0; jjj<addj; jjj++) {
00547                                                         for (int iii=0; iii<addi; iii++) {
00548                                                                 t+=image_data[l+iii+jjj*nx];
00549                                                         }
00550                                                 }
00551                                                 t/=addi*addi;
00552                                         }
00553                                         if (t <= rm) p = mingray;
00554                                         else if (t >= render_max) p = maxgray;
00555                                         else if (gamma!=1.0) {
00556                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00557                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00558 //                                              k = (int) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00559 //                                              k += mingray;
00560 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00561 //                                              k += mingray;
00562                                         }
00563                                         else {
00564                                                 p = (unsigned char) (gs * (t - render_min));
00565                                                 p += mingray;
00566                                         }
00567                                         data[i * asrgb + j * bpl] = p;
00568                                         if (hist) histd[p]++;
00569                                         l += addi;
00570                                         remx += addr;
00571                                         if (remx > scale_n) {
00572                                                 remx -= scale_n;
00573                                                 l++;
00574                                         }
00575                                 }
00576                                 l = br + addi * nx;
00577                                 remy += addr;
00578                                 if (remy > scale_n) {
00579                                         remy -= scale_n;
00580                                         l += nx;
00581                                 }
00582                         }
00583                 }
00584         }
00585 
00586         // this replicates r -> g,b
00587         if (asrgb==3 && !(flags&16)) {
00588                 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00589                         for (int i=xmin; i<xsize*3; i+=3) {
00590                                 data[i+j+1]=data[i+j+2]=data[i+j];
00591                         }
00592                 }
00593         }
00594         if (asrgb==4 && !(flags&16)) {
00595                 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00596                         for (int i=xmin; i<xsize*4; i+=4) {
00597                                 data[i+j+1]=data[i+j+2]=data[i+j+3]=data[i+j];
00598                                 data[i+j+3]=255;
00599                         }
00600                 }
00601         }
00602 
00603         EXITFUNC;
00604 
00605         // ok, ok, not the most efficient place to do this, but it works
00606         if (invy) {
00607                 int x,y;
00608                 char swp;
00609                 for (y=0; y<iysize/2; y++) {
00610                         for (x=0; x<ixsize; x++) {
00611                                 swp=ret[y*bpl+x];
00612                                 ret[y*bpl+x]=ret[(iysize-y-1)*bpl+x];
00613                                 ret[(iysize-y-1)*bpl+x]=swp;
00614                         }
00615                 }
00616         }
00617 
00618     //  return PyString_FromStringAndSize((const char*) data,iysize*bpl);
00619         if (flags&16) {
00620                 glDrawPixels(ixsize,iysize,GL_LUMINANCE,GL_UNSIGNED_BYTE,(const GLvoid *)ret.data());
00621         }
00622 
00623         return ret;
00624 }
00625 
00626 
00627 unsigned long GLUtil::get_isosurface_dl(MarchingCubes* mc, unsigned int tex_id,bool surface_face_z)
00628 {
00629         if ( mc->_isodl != 0 ) glDeleteLists(mc->_isodl,1);
00630 
00631 
00632         mc->calculate_surface();
00633         if ( surface_face_z ) mc->surface_face_z();
00634 #if MARCHING_CUBES_DEBUG
00635         cout << "There are " << ff.elem()/3 << " faces and " << pp.elem() << " points and " << nn.elem() << " normals to render in generate dl" << endl;
00636 #endif
00637         int maxf;
00638 //#ifdef        _WIN32
00639 //      glGetIntegerv(GL_MAX_ELEMENTS_INDICES_WIN,&maxf);
00640 //#else
00641         glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
00642 //#endif        //_WIN32
00643 #if MARCHING_CUBES_DEBUG
00644         int maxv;
00645         glGetIntegerv(GL_MAX_ELEMENTS_VERTICES,&maxv);
00646         cout << "Max vertices is " << maxv << " max indices is " << maxf << endl;
00647         cout << "Using OpenGL " << glGetString(GL_VERSION) << endl;
00648 #endif
00649 
00650         for (unsigned int i = 0; i < mc->ff.elem(); ++i ) mc->ff[i] /= 3;
00651 
00652         if ( maxf % 3 != 0 )
00653         {
00654                 maxf = maxf - (maxf%3);
00655         }
00656 
00657         if ( tex_id != 0 ) {
00658                 // Normalize the coordinates to be on the interval 0,1
00659                 mc->pp.mult3(1.0f/(float) mc->_emdata->get_xsize(), 1.0f/(float)mc->_emdata->get_ysize(), 1.0f/(float)mc->_emdata->get_zsize());
00660                 glEnableClientState(GL_TEXTURE_COORD_ARRAY);
00661                 glDisableClientState(GL_NORMAL_ARRAY);
00662                 glTexCoordPointer(3,GL_FLOAT,0,mc->pp.get_data());
00663         }
00664         else {
00665                 glEnableClientState(GL_NORMAL_ARRAY);
00666                 glDisableClientState(GL_TEXTURE_COORD_ARRAY);
00667                 glNormalPointer(GL_FLOAT,0,mc->nn.get_data());
00668         }
00669 
00670         glEnableClientState(GL_VERTEX_ARRAY);
00671         glVertexPointer(3,GL_FLOAT,0,mc->pp.get_data());
00672 
00673         mc->_isodl = glGenLists(1);
00674 
00675 #if MARCHING_CUBES_DEBUG
00676         int time0 = clock();
00677 #endif
00678         glNewList(mc->_isodl,GL_COMPILE);
00679 
00680         if ( tex_id != 0 ) {
00681                 glEnable(GL_TEXTURE_3D);
00682                 glBindTexture(GL_TEXTURE_3D, tex_id);
00683                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);
00684                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);
00685                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP);
00686                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
00687                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
00688                 glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
00689         }
00690         // Drawing range elements based on the output of glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
00691         // Saved about 60% of the time... drawRange should probably always be true
00692         bool drawRange = true;
00693         if ( drawRange == false ) {
00694                 glDrawElements(GL_TRIANGLES,mc->ff.elem(),GL_UNSIGNED_INT,mc->ff.get_data());
00695         } else {
00696                 for(unsigned int i = 0; i < mc->ff.elem(); i+=maxf)
00697                 {
00698                         if ( (i+maxf) > mc->ff.elem())
00699                                 glDrawElements(GL_TRIANGLES,mc->ff.elem()-i,GL_UNSIGNED_INT,&(mc->ff[i]));
00700                         else
00701                                 glDrawElements(GL_TRIANGLES,maxf,GL_UNSIGNED_INT,&(mc->ff[i]));
00702                                 // glDrawRangeElements is part of the extensions, we might want to experiment with its performance at some stage,
00703                                 // so please leave this code here, commented out. This is an either or situation, so if glDrawRangeElements is used,
00704                                 // glDrawElements above would have to be commented out.
00705                                 // glDrawRangeElements(GL_TRIANGLES,0,0,maxf,GL_UNSIGNED_INT,&ff[i]);
00706                 }
00707         }
00708 
00709         if ( tex_id != 0 ) glDisable(GL_TEXTURE_3D);
00710 
00711         glEndList();
00712 #if MARCHING_CUBES_DEBUG
00713         int time1 = clock();
00714         cout << "It took " << (time1-time0) << " " << (float)(time1-time0)/CLOCKS_PER_SEC << " to draw elements" << endl;
00715 #endif
00716         return mc->_isodl;
00717 }
00718 
00719 // This crap could be avoided and speed up(A lot) if we implemented an openGL shader........
00720 void GLUtil::glLoadMatrix(const Transform& xform)
00721 {
00722         vector<float> xformlist = xform.get_matrix_4x4();
00723         glLoadTransposeMatrixf(reinterpret_cast<GLfloat*>(&xformlist[0]));
00724 }
00725 
00726 void GLUtil::glMultMatrix(const Transform& xform)
00727 {
00728         vector<float> xformlist = xform.get_matrix_4x4();
00729         glMultTransposeMatrixf(reinterpret_cast<GLfloat*>(&xformlist[0]));
00730 }
00731 
00732 #endif // EMAN2_USING_OPENGL

Generated on Tue Jul 12 13:45:48 2011 for EMAN2 by  doxygen 1.4.7