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         if (emdata->is_complex()) {
00344                 if (dsx != -1) {
00345                         int l = y0 * nx;
00346                         for (int j = ymax; j >= ymin; j--) {
00347                                 int ll = x0;
00348                                 for (int i = xmin; i < xsize; i++) {
00349                                         if (l + ll > lmax || ll >= nx - 2) break;
00350 
00351                                         int k = 0;
00352                                         unsigned char p;
00353                                         if (ll >= nx / 2) {
00354                                                 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
00355                                                 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00356                                         }
00357                                         else k = nx * ny - (l + 2 * ll) - 2;
00358                                         if (k>=mid) k-=mid;             // These 2 lines handle the Fourier origin being in the corner, not the middle
00359                                         else k+=mid;
00360                                         float t = image_data[k];
00361                                         int ph;
00362                                         // in color mode
00363                                         if (flags&16 && asrgb>2) {
00364 //                                              if (l >= (ny - inv_scale) * nx) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;  // complex phase as integer 0-767;
00365                                                 if (ll >= nx / 2) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;        // complex phase as integer 0-767;
00366                                                 else ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;   // complex phase as integer 0-767;
00367                                         }
00368                                         if (t <= rm)  p = mingray;
00369                                         else if (t >= render_max) p = maxgray;
00370                                         else if (gamma!=1.0) {
00371                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00372                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00373 //                                              p = (unsigned char) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00374 //                                              p += mingray;
00375 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00376 //                                              k += mingray;
00377                                         }
00378                                         else {
00379                                                 p = (unsigned char) (gs * (t - render_min));
00380                                                 p += mingray;
00381                                         }
00382                                         // color rendering
00383                                         if (flags&16 && asrgb>2) {
00384                                                 if (ph<256) {
00385                                                         data[i * asrgb + j * bpl] = p*(255-ph)/256;
00386                                                         data[i * asrgb + j * bpl+1] = p*ph/256;
00387                                                         data[i * asrgb + j * bpl+2] = 0;
00388                                                 }
00389                                                 else if (ph<512) {
00390                                                         data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00391                                                         data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00392                                                         data[i * asrgb + j * bpl] = 0;
00393                                                 }
00394                                                 else {
00395                                                         data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00396                                                         data[i * asrgb + j * bpl] = p*(ph-512)/256;
00397                                                         data[i * asrgb + j * bpl+1] = 0;
00398                                                 }
00399                                         }
00400                                         else data[i * asrgb + j * bpl] = p;
00401                                         if (hist) histd[p]++;
00402                                         ll += dsx;
00403                                 }
00404                                 l += dsy;
00405                         }
00406                 }
00407                 else {
00408                         remy = 10;
00409                         int l = y0 * nx;
00410                         for (int j = ymax; j >= ymin; j--) {
00411                                 int br = l;
00412                                 remx = 10;
00413                                 int ll = x0;
00414                                 for (int i = xmin; i < xsize - 1; i++) {
00415                                         if (l + ll > lmax || ll >= nx - 2) {
00416                                                 break;
00417                                         }
00418                                         int k = 0;
00419                                         unsigned char p;
00420                                         if (ll >= nx / 2) {
00421                                                 if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
00422                                                 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00423                                         }
00424                                         else k = nx * ny - (l + 2 * ll) - 2;
00425                                         if (k>=mid) k-=mid;             // These 2 lines handle the Fourier origin being in the corner, not the middle
00426                                         else k+=mid;
00427 
00428                                         float t = image_data[k];
00429                                         // in color mode
00430                                         int ph;
00431                                         if (flags&16 && asrgb>2) {
00432                                                 if (l >= (ny * nx - nx)) ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;         // complex phase as integer 0-767;
00433                                                 else ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;   // complex phase as integer 0-767;
00434                                         }
00435                                         if (t <= rm)
00436                                                 p = mingray;
00437                                         else if (t >= render_max) {
00438                                                 p = maxgray;
00439                                         }
00440                                         else if (gamma!=1.0) {
00441                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00442                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00443 //                                              p = (unsigned char) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00444 //                                              p += mingray;
00445 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00446 //                                              k += mingray;
00447                                         }
00448                                         else {
00449                                                 p = (unsigned char) (gs * (t - render_min));
00450                                                 p += mingray;
00451                                         }
00452                                         if (flags&16 && asrgb>2) {
00453                                                 if (ph<256) {
00454                                                         data[i * asrgb + j * bpl] = p*(255-ph)/256;
00455                                                         data[i * asrgb + j * bpl+1] = p*ph/256;
00456                                                         data[i * asrgb + j * bpl+2] = 0;
00457                                                 }
00458                                                 else if (ph<512) {
00459                                                         data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00460                                                         data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00461                                                         data[i * asrgb + j * bpl] = 0;
00462                                                 }
00463                                                 else {
00464                                                         data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00465                                                         data[i * asrgb + j * bpl] = p*(ph-512)/256;
00466                                                         data[i * asrgb + j * bpl+1] = 0;
00467                                                 }
00468                                         }
00469                                         else data[i * asrgb + j * bpl] = p;
00470                                         if (hist) histd[p]++;
00471                                         ll += addi;
00472                                         remx += addr;
00473                                         if (remx > scale_n) {
00474                                                 remx -= scale_n;
00475                                                 ll++;
00476                                         }
00477                                 }
00478                                 l = br + addi * nx;
00479                                 remy += addr;
00480                                 if (remy > scale_n) {
00481                                         remy -= scale_n;
00482                                         l += nx;
00483                                 }
00484                         }
00485                 }
00486         }
00487         else {
00488                 if (dsx != -1) {
00489                         int l = x0 + y0 * nx;
00490                         for (int j = ymax; j >= ymin; j--) {
00491                                 int br = l;
00492                                 for (int i = xmin; i < xsize; i++) {
00493                                         if (l > lmax) {
00494                                                 break;
00495                                         }
00496                                         int k = 0;
00497                                         unsigned char p;
00498                                         float t;
00499                                         if (dsx==1) t=image_data[l];
00500                                         else {                                          // This block does local pixel averaging for nicer reduced views
00501                                                 t=0;
00502                                                 for (int iii=0; iii<dsx; iii++) {
00503                                                         for (int jjj=0; jjj<dsy; jjj+=nx) {
00504                                                                 t+=image_data[l+iii+jjj];
00505                                                         }
00506                                                 }
00507                                                 t/=dsx*(dsy/nx);
00508                                         }
00509 
00510                                         if (t <= rm) p = mingray;
00511                                         else if (t >= render_max) p = maxgray;
00512                                         else if (gamma!=1.0) {
00513                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00514                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00515 //                                              k = (int) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00516 //                                              k += mingray;
00517 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00518 //                                              k += mingray;
00519                                         }
00520                                         else {
00521                                                 p = (unsigned char) (gs * (t - render_min));
00522                                                 p += mingray;
00523                                         }
00524                                         data[i * asrgb + j * bpl] = p;
00525                                         if (hist) histd[p]++;
00526                                         l += dsx;
00527                                 }
00528                                 l = br + dsy;
00529                         }
00530                 }
00531                 else {
00532                         remy = 10;
00533                         int l = x0 + y0 * nx;
00534                         for (int j = ymax; j >= ymin; j--) {
00535                                 int addj = addi;
00536                                 // There seems to be some overflow issue happening
00537                                 // where the statement if (l > lmax) break (below) doesn't work
00538                                 // because the loop that iterates jjj can inadvertantly go out of bounds
00539                                 if (( l + addi*nx ) >= nxy ) {
00540                                         addj = (nxy-l)/nx;
00541                                         if (addj <= 0) continue;
00542                                 }
00543                                 int br = l;
00544                                 remx = 10;
00545                                 for (int i = xmin; i < xsize; i++) {
00546                                         if (l > lmax) break;
00547                                         int k = 0;
00548                                         unsigned char p;
00549                                         float t;
00550                                         if (addi<=1) t = image_data[l];
00551                                         else {                                          // This block does local pixel averaging for nicer reduced views
00552                                                 t=0;
00553                                                 for (int jjj=0; jjj<addj; jjj++) {
00554                                                         for (int iii=0; iii<addi; iii++) {
00555                                                                 t+=image_data[l+iii+jjj*nx];
00556                                                         }
00557                                                 }
00558                                                 t/=addi*addi;
00559                                         }
00560                                         if (t <= rm) p = mingray;
00561                                         else if (t >= render_max) p = maxgray;
00562                                         else if (gamma!=1.0) {
00563                                                 k=(int)(gs2 * (t-render_min));          // map float value to 0-4096 range
00564                                                 p = gammamap[k];                                        // apply gamma using precomputed gamma map
00565 //                                              k = (int) (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma);
00566 //                                              k += mingray;
00567 //                                              k = static_cast<int>( (maxgray-mingray)*pow((gs2 * (t - render_min)),gamma) );
00568 //                                              k += mingray;
00569                                         }
00570                                         else {
00571                                                 p = (unsigned char) (gs * (t - render_min));
00572                                                 p += mingray;
00573                                         }
00574                                         data[i * asrgb + j * bpl] = p;
00575                                         if (hist) histd[p]++;
00576                                         l += addi;
00577                                         remx += addr;
00578                                         if (remx > scale_n) {
00579                                                 remx -= scale_n;
00580                                                 l++;
00581                                         }
00582                                 }
00583                                 l = br + addi * nx;
00584                                 remy += addr;
00585                                 if (remy > scale_n) {
00586                                         remy -= scale_n;
00587                                         l += nx;
00588                                 }
00589                         }
00590                 }
00591         }
00592 
00593         // this replicates r -> g,b
00594         if (asrgb==3 && !(flags&16)) {
00595                 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00596                         for (int i=xmin; i<xsize*3; i+=3) {
00597                                 data[i+j+1]=data[i+j+2]=data[i+j];
00598                         }
00599                 }
00600         }
00601         if (asrgb==4 && !(flags&16)) {
00602                 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00603                         for (int i=xmin; i<xsize*4; i+=4) {
00604                                 data[i+j+1]=data[i+j+2]=data[i+j+3]=data[i+j];
00605                                 data[i+j+3]=255;
00606                         }
00607                 }
00608         }
00609 
00610         EXITFUNC;
00611 
00612         // ok, ok, not the most efficient place to do this, but it works
00613         if (invy) {
00614                 int x,y;
00615                 char swp;
00616                 for (y=0; y<iysize/2; y++) {
00617                         for (x=0; x<ixsize; x++) {
00618                                 swp=ret[y*bpl+x];
00619                                 ret[y*bpl+x]=ret[(iysize-y-1)*bpl+x];
00620                                 ret[(iysize-y-1)*bpl+x]=swp;
00621                         }
00622                 }
00623         }
00624 
00625     //  return PyString_FromStringAndSize((const char*) data,iysize*bpl);
00626         if (flags&16) {
00627                 glDrawPixels(ixsize,iysize,GL_LUMINANCE,GL_UNSIGNED_BYTE,(const GLvoid *)ret.data());
00628         }
00629 
00630         return ret;
00631 }
00632 
00633 //DEPRECATED
00634 unsigned long GLUtil::get_isosurface_dl(MarchingCubes* mc, unsigned int tex_id,bool surface_face_z, bool recontour)
00635 {
00636         if ( mc->_isodl != 0 ) glDeleteLists(mc->_isodl,1);
00637                 
00638         if ( recontour ) mc->calculate_surface();
00639         if ( surface_face_z ) mc->surface_face_z();
00640         if( mc->getRGBmode() ) mc->color_vertices();
00641         
00642 #if MARCHING_CUBES_DEBUG
00643         cout << "There are " << ff.elem()/3 << " faces and " << pp.elem() << " points and " << nn.elem() << " normals to render in generate dl" << endl;
00644 #endif
00645         int maxf;
00646 //#ifdef        _WIN32
00647 //      glGetIntegerv(GL_MAX_ELEMENTS_INDICES_WIN,&maxf);
00648 //#else
00649         glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
00650 //#endif        //_WIN32
00651 #if MARCHING_CUBES_DEBUG
00652         int maxv;
00653         glGetIntegerv(GL_MAX_ELEMENTS_VERTICES,&maxv);
00654         cout << "Max vertices is " << maxv << " max indices is " << maxf << endl;
00655         cout << "Using OpenGL " << glGetString(GL_VERSION) << endl;
00656 #endif
00657 
00658         if ( recontour ) for (unsigned int i = 0; i < mc->ff.elem(); ++i ) mc->ff[i] /= 3;
00659 
00660         if ( maxf % 3 != 0 )
00661         {
00662                 maxf = maxf - (maxf%3);
00663         }
00664 
00665         if ( tex_id != 0 ) {
00666                 // Normalize the coordinates to be on the interval 0,1
00667                 mc->pp.mult3(1.0f/(float) mc->_emdata->get_xsize(), 1.0f/(float)mc->_emdata->get_ysize(), 1.0f/(float)mc->_emdata->get_zsize());
00668                 glEnableClientState(GL_TEXTURE_COORD_ARRAY);
00669                 glDisableClientState(GL_NORMAL_ARRAY);
00670                 glTexCoordPointer(3,GL_FLOAT,0,mc->pp.get_data());
00671         }
00672         else {
00673                 glEnableClientState(GL_NORMAL_ARRAY);
00674                 glDisableClientState(GL_TEXTURE_COORD_ARRAY);
00675                 glNormalPointer(GL_FLOAT,0, mc->nn.get_data());
00676         }
00677 
00678         glEnableClientState(GL_VERTEX_ARRAY);
00679         glVertexPointer(3,GL_FLOAT,0, mc->pp.get_data());
00680         
00681         if (mc->getRGBmode()){
00682                 glEnableClientState(GL_COLOR_ARRAY);
00683                 glColorPointer(3,GL_FLOAT,0, mc->cc.get_data());
00684         }
00685         
00686         mc->_isodl = glGenLists(1);
00687         
00688 #if MARCHING_CUBES_DEBUG
00689         int time0 = clock();
00690 #endif
00691         
00692         glNewList(mc->_isodl,GL_COMPILE);
00693 
00694         if ( tex_id != 0 ) {
00695                 glEnable(GL_TEXTURE_3D);
00696                 glBindTexture(GL_TEXTURE_3D, tex_id);
00697                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);
00698                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);
00699                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP);
00700                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
00701                 glTexParameterf(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
00702                 glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
00703         }
00704         
00705         // Drawing range elements based on the output of glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
00706         // Saved about 60% of the time... drawRange should probably always be true
00707         bool drawRange = true;
00708         if ( drawRange == false ) {
00709                 glDrawElements(GL_TRIANGLES,mc->ff.elem(),GL_UNSIGNED_INT, mc->ff.get_data());
00710         } else {
00711                 for(unsigned int i = 0; i < mc->ff.elem(); i+=maxf)
00712                 {
00713                         if ( (i+maxf) > mc->ff.elem())
00714                                 glDrawElements(GL_TRIANGLES,mc->ff.elem()-i,GL_UNSIGNED_INT,&(mc->ff[i]));
00715                         else
00716                                 glDrawElements(GL_TRIANGLES,maxf,GL_UNSIGNED_INT,&(mc->ff[i]));
00717                                 // glDrawRangeElements is part of the extensions, we might want to experiment with its performance at some stage,
00718                                 // so please leave this code here, commented out. This is an either or situation, so if glDrawRangeElements is used,
00719                                 // glDrawElements above would have to be commented out.
00720                                 // glDrawRangeElements(GL_TRIANGLES,0,0,maxf,GL_UNSIGNED_INT,&ff[i]);
00721                 }
00722         }
00723 
00724         if ( tex_id != 0 ) glDisable(GL_TEXTURE_3D);
00725 
00726         glEndList();
00727         
00728 #if MARCHING_CUBES_DEBUG
00729         int time1 = clock();
00730         cout << "It took " << (time1-time0) << " " << (float)(time1-time0)/CLOCKS_PER_SEC << " to draw elements" << endl;
00731 #endif
00732         return mc->_isodl;
00733 }
00734 
00735 void GLUtil::contour_isosurface(MarchingCubes* mc)
00736 {
00737         mc->calculate_surface();
00738         //What does this do???
00739         for (unsigned int i = 0; i < mc->ff.elem(); ++i ) mc->ff[i] /= 3;
00740         //Need to rebind data (to the GPU)
00741         mc->needtobind = true;
00742 }
00743 
00744 void GLUtil::render_using_VBOs(MarchingCubes* mc, unsigned int tex_id,bool surface_face_z)
00745 {
00746         // In current version Texture is not supported b/c it is not used... EVER
00747         
00748 #ifdef _WIN32
00749         typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
00750         PFNGLGENBUFFERSPROC glGenBuffers;
00751         glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
00752 
00753         typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
00754         PFNGLISBUFFERPROC glIsBuffer;
00755         glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
00756 #endif  //_WIN32
00757 
00758         if ( surface_face_z ) mc->surface_face_z();
00759         
00760         //Bug in OpenGL, sometimes glGenBuffers doesn't work. Try again, if we still fail then return...
00761         if (!glIsBuffer(mc->buffer[0])) glGenBuffers(4, mc->buffer);
00762 
00763         //whenever something changes, like color mode or color scale (or threshold), we need to recolor
00764         if( mc->getRGBmode() && (mc->rgbgenerator.getNeedToRecolor() || mc->needtobind)){
00765                 mc->color_vertices();
00766                 mc->needtobind = true;
00767         }
00768 
00769         int maxf;
00770         glGetIntegerv(GL_MAX_ELEMENTS_INDICES,&maxf);
00771 
00772         if ( maxf % 3 != 0 )
00773         {
00774                 maxf = maxf - (maxf%3);
00775         }
00776 
00777 #ifdef _WIN32
00778         typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
00779         PFNGLBINDBUFFERPROC glBindBuffer;
00780         glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
00781 
00782         typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size, const GLvoid *data, GLenum usage);
00783         PFNGLBUFFERDATAPROC glBufferData;
00784         glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
00785 #endif  //_WIN32
00786 
00787         //Normal vectors
00788         glEnableClientState(GL_NORMAL_ARRAY);
00789         glDisableClientState(GL_TEXTURE_COORD_ARRAY);
00790         glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[2]);
00791         if (mc->needtobind) glBufferData(GL_ARRAY_BUFFER, 4*mc->nn.elem(), mc->nn.get_data(), GL_STATIC_DRAW);
00792         glNormalPointer(GL_FLOAT,0,0);
00793 
00794         //Vertex vectors
00795         glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[0]);
00796         if (mc->needtobind) glBufferData(GL_ARRAY_BUFFER, 4*mc->pp.elem(), mc->pp.get_data(), GL_STATIC_DRAW);
00797         glEnableClientState(GL_VERTEX_ARRAY);
00798         glVertexPointer(3,GL_FLOAT,0,0);
00799 
00800         //Color vectors
00801         if (mc->getRGBmode()){
00802                 glEnableClientState(GL_COLOR_ARRAY);
00803                 glBindBuffer(GL_ARRAY_BUFFER, mc->buffer[3]);
00804                 if (mc->needtobind) glBufferData(GL_ARRAY_BUFFER, 4*mc->cc.elem(), mc->cc.get_data(), GL_STATIC_DRAW);
00805                 glColorPointer(3,GL_FLOAT,0, 0);
00806         }
00807         
00808         //Indices
00809         glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, mc->buffer[1]);
00810         if (mc->needtobind) glBufferData(GL_ELEMENT_ARRAY_BUFFER, 4*mc->ff.elem(), mc->ff.get_data(), GL_STATIC_DRAW);
00811         
00812         // This lets us know if buffers an not implmemted
00813         if (!glIsBuffer(mc->buffer[0])){
00814                 cout << "Can't Generate GL Vertex Buffers. glGenBuffer error" << endl;
00815                 return;
00816         }
00817         
00818         //finally draw the elemenets
00819         glDrawElements(GL_TRIANGLES,mc->ff.elem(),GL_UNSIGNED_INT,0);
00820         // No longer need to bind data to the GPU
00821         mc->needtobind = false;
00822 
00823 }
00824 
00825 // This crap could be avoided and speed up(A lot) if we implemented an openGL shader........
00826 void GLUtil::glLoadMatrix(const Transform& xform)
00827 {
00828 #ifdef _WIN32
00829         typedef void (APIENTRYP PFNGLLOADTRANSPOSEMATRIXFPROC) (const GLfloat *m);
00830         PFNGLLOADTRANSPOSEMATRIXFPROC glLoadTransposeMatrixf = NULL;
00831         glLoadTransposeMatrixf = (PFNGLLOADTRANSPOSEMATRIXFPROC) wglGetProcAddress("glLoadTransposeMatrixf");
00832 #endif  //_WIN32
00833 
00834         vector<float> xformlist = xform.get_matrix_4x4();
00835         glLoadTransposeMatrixf(reinterpret_cast<GLfloat*>(&xformlist[0]));
00836 }
00837 
00838 void GLUtil::glMultMatrix(const Transform& xform)
00839 {
00840 #ifdef _WIN32
00841         typedef void (APIENTRYP PFNGLMULTTRANSPOSEMATRIXFPROC) (const GLfloat *m);
00842         PFNGLMULTTRANSPOSEMATRIXFPROC glMultTransposeMatrixf = NULL;
00843         glMultTransposeMatrixf = (PFNGLMULTTRANSPOSEMATRIXFPROC) wglGetProcAddress("glMultTransposeMatrixf");
00844 #endif  //_WIN32
00845 
00846         vector<float> xformlist = xform.get_matrix_4x4();
00847         glMultTransposeMatrixf(reinterpret_cast<GLfloat*>(&xformlist[0]));
00848 }
00849 
00850 //This draw a bounding box, or any box 
00851 void GLUtil::glDrawBoundingBox(float width, float height, float depth)
00852 {
00853         
00854         float w2 = width/2.0f;
00855         float h2 = height/2.0f;
00856         float d2 = depth/2.0f;
00857         
00858         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};
00859         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};
00860         
00861 #ifdef _WIN32
00862         typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
00863         PFNGLGENBUFFERSPROC glGenBuffers;
00864         glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
00865 
00866         typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
00867         PFNGLISBUFFERPROC glIsBuffer;
00868         glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
00869 
00870         typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
00871         PFNGLBINDBUFFERPROC glBindBuffer;
00872         glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
00873 
00874         typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size, const GLvoid *data, GLenum usage);
00875         PFNGLBUFFERDATAPROC glBufferData;
00876         glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
00877 #endif  //_WIN32
00878 
00879         if ( glIsBuffer(GLUtil::buffer[0])  == GL_FALSE ){
00880                 glGenBuffers(2, GLUtil::buffer);
00881         }
00882         
00883         // Could use dirty bit here but not worth my time to implment
00884         glBindBuffer(GL_ARRAY_BUFFER, GLUtil::buffer[0]);
00885         glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);
00886         glEnableClientState(GL_VERTEX_ARRAY);
00887         glVertexPointer(3,GL_FLOAT,0,0);
00888         
00889         glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, GLUtil::buffer[1]);
00890         glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(indices), indices, GL_STATIC_DRAW);
00891         glDrawElements(GL_LINES,24,GL_UNSIGNED_INT,0);
00892         
00893 }
00894 
00895 void GLUtil::glDrawDisk(float radius, int spokes)
00896 {
00897 #ifdef _WIN32
00898         typedef void (APIENTRYP PFNGLGENBUFFERSPROC) (GLsizei n, GLuint *buffers);
00899         PFNGLGENBUFFERSPROC glGenBuffers;
00900         glGenBuffers = (PFNGLGENBUFFERSPROC) wglGetProcAddress("glGenBuffers");
00901 
00902         typedef GLboolean (APIENTRYP PFNGLISBUFFERPROC) (GLuint buffer);
00903         PFNGLISBUFFERPROC glIsBuffer;
00904         glIsBuffer = (PFNGLISBUFFERPROC) wglGetProcAddress("glIsBuffer");
00905 
00906         typedef void (APIENTRYP PFNGLBINDBUFFERPROC) (GLenum target, GLuint buffer);
00907         PFNGLBINDBUFFERPROC glBindBuffer;
00908         glBindBuffer = (PFNGLBINDBUFFERPROC) wglGetProcAddress("glBindBuffer");
00909 
00910         typedef void (APIENTRYP PFNGLBUFFERDATAPROC) (GLenum target, GLsizeiptr size, const GLvoid *data, GLenum usage);
00911         PFNGLBUFFERDATAPROC glBufferData;
00912         glBufferData = (PFNGLBUFFERDATAPROC) wglGetProcAddress("glBufferData");
00913 #endif  //_WIN32
00914 
00915         //This code is experimental
00916         int arraysize = 4*pow((double)2,(double)spokes);
00917         int sideofarray = pow((double)2,(double)spokes);
00918         
00919 //      float vertices[3*arraysize + 3];
00920         vector<float> vertices(3*arraysize + 3);
00921         
00922         //last vertex is center
00923         vertices[3*arraysize] = 0.0;
00924         vertices[3*arraysize+1] = 0.0;
00925         vertices[3*arraysize+2] = 0.0;
00926         
00927         vertices[0] = 0.0;
00928         vertices[1] = radius;
00929         vertices[2] = 0.0;
00930                 
00931         vertices[sideofarray*3] = radius;
00932         vertices[sideofarray*3 + 1] = 0.0;
00933         vertices[sideofarray*3 + 2] = 0.0;
00934         
00935         vertices[sideofarray*6] = 0.0;
00936         vertices[sideofarray*6 + 1] = -radius;
00937         vertices[sideofarray*6 + 2] = 0.0;
00938         
00939         vertices[sideofarray*9] = -radius;
00940         vertices[sideofarray*9 + 1] = 0.0;
00941         vertices[sideofarray*9 + 2] = 0.0;
00942         
00943         //This could aslo be implemented recusivly
00944         for(int step = 0; step < spokes; step++){
00945                 //starting location
00946                 int x = sideofarray/pow(2.0,(double)(step+1));
00947                 for(int i = 1; i <= 4*pow(2.0,(double)step); i++ ){
00948                         
00949                         // take the necessary steps
00950                         int index =  x + 2*x*(i-1);
00951                         int index_f = (index + x) % arraysize;
00952                         int index_i = index - x;
00953                         cout << index << " " << index_f << " " << index_i << endl;
00954                         
00955                         //need to resclae length to that of radius
00956                         vertices[index_f*3] = (vertices[index_f*3] - vertices[index_i*3])/2.0f;
00957                         vertices[index_f*3 + 1] = (vertices[index_f*3 + 1] - vertices[index_i*3 + 1])/2.0f;
00958                         vertices[index_f*3 + 2] = (vertices[index_f*3 + 2] - vertices[index_i*3 + 2])/2.0f;
00959                 }
00960         }
00961         
00962         //GL stuff
00963         if ( glIsBuffer(GLUtil::buffer[0])  == GL_FALSE ){
00964                 glGenBuffers(2, GLUtil::buffer);
00965         }
00966         
00967         // Could use dirty bit here but not worth my time to implment
00968         glBindBuffer(GL_ARRAY_BUFFER, GLUtil::buffer[0]);
00969         //glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);
00970         glBufferData(GL_ARRAY_BUFFER, vertices.size(), &(vertices[0]), GL_STATIC_DRAW);
00971         glEnableClientState(GL_VERTEX_ARRAY);
00972         glVertexPointer(3,GL_FLOAT,0,0);
00973         
00974         //Code to select indices
00975         
00976         
00977 }
00978 #endif // EMAN2_USING_OPENGL

Generated on Thu May 3 10:06:24 2012 for EMAN2 by  doxygen 1.4.7