Changeset 378 for trunk/gui/vizservers

Timestamp:

Mar 27, 2006, 6:56:57 AM (18 years ago)

Author:

mmc

Message:

Fixed loading of SQUALID-2D data.

Location:

trunk/gui/vizservers/nanovis

Files:

• ParticleSystem.cpp (modified) (6 diffs)
• ParticleSystem.h (modified) (3 diffs)
• nanovis.cpp (modified) (21 diffs)
• nanovis.h (modified) (2 diffs)
• shaders/one_volume.cg (modified) (1 diff)
• shaders/render_vel.cg (modified) (1 diff)
• shaders/update_pos.cg (modified) (2 diffs)

trunk/gui/vizservers/nanovis/ParticleSystem.cpp

@@ -22 +22 @@
 
 
-ParticleSystem::ParticleSystem(int w, int h, CGcontext context, NVISid volume){
+ParticleSystem::ParticleSystem(int w, int h, CGcontext context, NVISid volume, float scale_x, 
+                float scale_y, float scale_z){
+
+  fprintf(stderr, "%f, %f, %f\n", scale_x, scale_y, scale_z);
 
   psys_width = w;
@@ -30 +33 @@
   reborn = true;
   flip = true;
-  max_life = 30;
+  max_life = 100;
 
   data = (Particle*) malloc(w*h*sizeof(Particle));
@@ -66 +69 @@
   m_vel_tex_param = cgGetNamedParameter(m_pos_fprog, "vel_tex");
   m_pos_tex_param = cgGetNamedParameter(m_pos_fprog, "pos_tex");
+  m_scale_param = cgGetNamedParameter(m_pos_fprog, "scale");
   cgGLSetTextureParameter(m_vel_tex_param, volume);
+  cgGLSetParameter3f(m_scale_param, scale_x, scale_y, scale_z);
 
   fprintf(stderr, "init_psys\n");
@@ -131 +136 @@
 
      cgGLBindProgram(m_pos_fprog);
-     cgGLSetParameter1f(m_pos_timestep_param, 0.1);
+     cgGLSetParameter1f(m_pos_timestep_param, 0.05);
      cgGLEnableTextureParameter(m_vel_tex_param);
      cgGLSetTextureParameter(m_pos_tex_param, psys_tex[0]);
@@ -159 +164 @@
 
      cgGLBindProgram(m_pos_fprog);
-     cgGLSetParameter1f(m_pos_timestep_param, 0.1);
+     cgGLSetParameter1f(m_pos_timestep_param, 0.05);
      cgGLEnableTextureParameter(m_vel_tex_param);
      cgGLSetTextureParameter(m_pos_tex_param, psys_tex[1]);
@@ -200 +205 @@
   glDisable(GL_BLEND);
 
-  glPointSize(0.5);
-  glColor4f(.8,.8,.8,1.);
+  //glPointSize(0.5);
+  glPointSize(1.5);
+  glColor4f(.1,.0,.3,1.);
 
   m_vertex_array->SetPointer(0);

trunk/gui/vizservers/nanovis/ParticleSystem.h

@@ -34 +34 @@
 
   Particle(){};
-  Particle(float _x, float _y, float _z, float _life):
+  Particle(float _x, float _y, float _z, float _life) :
    x(_x), y(_y), z(_z), aux(_life){}
 };
@@ -56 +56 @@
   CGcontext m_g_context;
   CGprogram m_pos_fprog;
-  CGparameter m_vel_tex_param, m_pos_tex_param;
+  CGparameter m_vel_tex_param, m_pos_tex_param, m_scale_param;
   CGparameter m_pos_timestep_param, m_pos_spherePos_param;
 
@@ -63 +63 @@
   int psys_height;
 
-  ParticleSystem(int w, int h, CGcontext context, NVISid vector_field);
+  ParticleSystem(int w, int h, CGcontext context, NVISid vector_field, 
+                  float scalex, float scaley, float scalez);
   ~ParticleSystem();
   void initialize(Particle* data);

trunk/gui/vizservers/nanovis/nanovis.cpp

@@ -23 +23 @@
 #include "RpFieldPrism3D.h"
 
+// forward declarations
+void init_particles();
+void makePatterns();
+void get_slice_vectors();
+
 ParticleSystem* psys;
 float psys_x=0.4, psys_y=0, psys_z=0;
@@ -116 +121 @@
 void load_volume(int index, int width, int height, int depth, int n_component, float* data);
 
-
-//load a vector field from file
-void init_vector_field(){
-
-  FILE* fp = fopen("./data/tornado_64x64x64.raw", "rb");
-  assert(fp!=0);
-
-  float* data = new float[64*64*64*3];
-  for(int j=0; j<64*64*64*3; j++){
-       fread(data + j, 4, 1, fp);
-       //fprintf(stderr, "%f ", data[3*j+i]);
-  }
-  fclose(fp);
-
-  //normalize data
-  float _max = FLT_MIN;
-  float _min = FLT_MAX;
-
-  for(int j=0; j<64*64*64; j++){
-    float mag = sqrt(data[j*3]*data[j*3] + data[j*3+1]*data[j*3+1] + data[j*3+2]*data[j*3+2]);
-    if(mag > _max)
-      _max = mag;
-    if(mag < _min) 
-      _min = mag;
-  }
-
-  fprintf(stderr, "max=%f, min=%f\n", _max, _min);
-
-  float scale = _max - _min;
-  
-  for(int j=0; j<64*64*64*3; j++){
-    data[j]=data[j]/scale;
+/* Load a 3D vector field from a dx-format file
+ */
+void
+load_vector_file(int index, char *fname) {
+    int dummy, nx, ny, nz, nxy, npts;
+    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
+    char line[128], type[128], *start;
+    std::ifstream fin(fname);
+
+    do {
+        fin.getline(line,sizeof(line)-1);
+        for (start=&line[0]; *start == ' ' || *start == '\t'; start++)
+            ;  // skip leading blanks
+
+        if (*start != '#') {  // skip comment lines
+            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
+                // found grid size
+            }
+            else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
+                // found origin
+            }
+            else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
+                // found one of the delta lines
+                if (ddx != 0.0) { dx = ddx; }
+                else if (ddy != 0.0) { dy = ddy; }
+                else if (ddz != 0.0) { dz = ddz; }
+            }
+            else if (sscanf(start, "object %d class array type %s shape 3 rank 1 items %d data follows", &dummy, type, &npts) == 3) {
+                if (npts != nx*ny*nz) {
+                    std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << npts << " points" << std::endl;
+                    return;
+                }
+                break;
+            }
+        }
+    } while (!fin.eof());
+
+    // read data points
+    if (!fin.eof()) {
+        Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
+        Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
+        Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
+        Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
+        Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
+        Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);
+
+        double vx, vy, vz;
+        int nread = 0;
+        for (int ix=0; ix < nx; ix++) {
+            for (int iy=0; iy < ny; iy++) {
+                for (int iz=0; iz < nz; iz++) {
+                    if (fin.eof() || nread > npts) {
+                        break;
+                    }
+                    fin.getline(line,sizeof(line)-1);
+                    if (sscanf(line, "%lg %lg %lg", &vx, &vy, &vz) == 3) {
+                        int nindex = iz*nx*ny + iy*nx + ix;
+                        xfield.define(nindex, vx);
+                        yfield.define(nindex, vy);
+                        zfield.define(nindex, vz);
+                        nread++;
+                    }
+                }
+            }
+        }
+
+        // make sure that we read all of the expected points
+        if (nread != nx*ny*nz) {
+            std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << nread << " points" << std::endl;
+            return;
+        }
+
+        // figure out a good mesh spacing
+        int nsample = 30;
+        dx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
+        dy = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
+        dz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);
+        double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
+
+        nx = (int)ceil(dx/dmin);
+        ny = (int)ceil(dy/dmin);
+        nz = (int)ceil(dz/dmin);
 
 #ifndef NV40
-    //cut negative values
-    if(data[j]<0)
-      data[j]=0;
+        nx = pow(2.0, ceil(log10((double)nx)/log10(2.0)));  // must be an even power of 2
+        ny = pow(2.0, ceil(log10((double)ny)/log10(2.0)));
+        nz = pow(2.0, ceil(log10((double)nz)/log10(2.0)));
 #endif
-  }
-
-  /*
-  for(int j=0; j<64*64*64; j++)
-    fprintf(stderr, "[%f, %f, %f]\t", data[j*3], data[j*3+1], data[j*3+2]);
-  */
-
-  //load the vector field as a volume
-  load_volume(0, 64, 64, 64, 3, data);
-  n_volumes++;
-  delete[] data;
-
-  fprintf(stderr, "init_vector_field\n");
-}
+
+        float *data = new float[3*nx*ny*nz];
+
+        std::cout << "generating " << nx << "x" << ny << "x" << nz << " = " << nx*ny*nz << " points" << std::endl;
+
+        // generate the uniformly sampled data that we need for a volume
+        double vmin = 0.0;
+        double vmax = 0.0;
+        int ngen = 0;
+        for (int iz=0; iz < nz; iz++) {
+            double zval = z0 + iz*dmin;
+            for (int iy=0; iy < ny; iy++) {
+                double yval = y0 + iy*dmin;
+                for (int ix=0; ix < nx; ix++) {
+                    double xval = x0 + ix*dmin;
+
+                    vx = xfield.value(xval,yval,zval);
+                    vy = yfield.value(xval,yval,zval);
+                    vz = zfield.value(xval,yval,zval);
+                    //vx = 1;
+                    //vy = 1;
+                    vz = 0;
+                    double vm = sqrt(vx*vx + vy*vy + vz*vz);
+
+                    if (vm < vmin) { vmin = vm; }
+                    if (vm > vmax) { vmax = vm; }
+
+                    data[ngen++] = vx;
+                    data[ngen++] = vy;
+                    data[ngen++] = vz;
+                }
+            }
+        }
+
+        ngen = 0;
+        for (ngen=0; ngen < npts; ngen++) {
+            data[ngen] = (data[ngen]/(2.0*vmax) + 0.5);
+        }
+
+        load_volume(index, nx, ny, nz, 3, data);
+        delete [] data;
+    } else {
+        std::cerr << "WARNING: data not found in file " << fname << std::endl;
+    }
+}
+
 
 /* Load a 3D volume from a dx-format file
@@ -270 +357 @@
             double dval;
             int nread = 0;
-            while (!fin.eof() && nread < npts) {
-                if (!(fin >> dval).fail()) {
-                    field.define(nread++, dval);
+            for (int ix=0; ix < nx; ix++) {
+                for (int iy=0; iy < ny; iy++) {
+                    for (int iz=0; iz < nz; iz++) {
+                        if (fin.eof() || nread > npts) {
+                            break;
+                        }
+                        fin.getline(line,sizeof(line)-1);
+                        if (sscanf(line, "%lg", &dval) == 1) {
+                            int nindex = iz*nx*ny + iy*nx + ix;
+                            field.define(nindex, dval);
+                            nread++;
+                        }
+                    }
                 }
             }
@@ -532 +629 @@
     m_mvi_vert_std_param = cgGetNamedParameter(m_vert_std_vprog, "modelViewInv");
 
+/*
     m_pos_fprog = loadProgram(g_context, CG_PROFILE_FP30, CG_SOURCE, "./shaders/update_pos.cg");
     m_pos_timestep_param  = cgGetNamedParameter(m_pos_fprog, "timestep");
@@ -537 +635 @@
     m_pos_tex_param = cgGetNamedParameter(m_pos_fprog, "pos_tex");
     cgGLSetTextureParameter(m_vel_tex_param, volume[0]->id);
+*/
+  
 
     m_render_vel_fprog = loadProgram(g_context, CG_PROFILE_FP30, CG_SOURCE, "./shaders/render_vel.cg");
@@ -602 +702 @@
       int index = i + psys->psys_height*j;
       bool particle = rand() % 256 > 100; 
+      particle = true;
       if(particle)
       {
-        data[4*index] = psys_x;
-        data[4*index+1]= i/float(psys->psys_width);
-        data[4*index+2]= j/float(psys->psys_height);
-        data[4*index+3]= 30;    
+        //assign any location (x,y,z) in range [0,1]
+        data[4*index] = 0.5;
+        data[4*index+1]= j/float(psys->psys_height);
+        data[4*index+2]= lic_slice_z; //lic_slice_z; //i/float(psys->psys_width);
+        data[4*index+3]= 30; //shorter life span, quicker iterations    
       }
       else
@@ -627 +729 @@
 
 
-void makePatterns();
-void get_slice_vectors();
 
 
@@ -680 +780 @@
    //load_volume_file(0, "./data/test2.dx");
 
+   load_vector_file(0, "./data/J-wire-vec.dx");
+   load_volume_file(1, "./data/mu-wire-3d.dx");
+
    init_fbo();  //frame buffer objects
    init_cg();   //init cg shaders
@@ -686 +789 @@
    switch_shader(cur_shader);
 
-   //psys = new ParticleSystem(NMESH, NMESH, g_context, volume[0]->id);
-   //init_particles();  //fill initial particles
-
-   //get_slice_vectors();
+   /*
+   psys = new ParticleSystem(NMESH, NMESH, g_context, volume[0]->id,
+                   1./volume[0]->aspect_ratio_width,
+                   1./volume[0]->aspect_ratio_height,
+                   1./volume[0]->aspect_ratio_depth);
+   */
+   psys = new ParticleSystem(NMESH, NMESH, g_context, volume[0]->id, 0.25, 1., 1.);
+
+   init_particles();    //fill initial particles
+
+   get_slice_vectors();
 }
 
@@ -1085 +1195 @@
    cgGLEnableProfile(CG_PROFILE_FP30);
    glBegin(GL_QUADS);
-    glTexCoord3f(0., 0., slice_z); glVertex2f(0.,         0.);
-    glTexCoord3f(1., 0., slice_z); glVertex2f((float)NMESH,     0.);
-    glTexCoord3f(1., 1., slice_z); glVertex2f((float)NMESH, (float)NMESH);
-    glTexCoord3f(0., 1., slice_z); glVertex2f((float)0.,    (float)NMESH);
+    glTexCoord3f(0., 0., lic_slice_z); glVertex2f(0.,         0.);
+    glTexCoord3f(1., 0., lic_slice_z); glVertex2f((float)NMESH,     0.);
+    glTexCoord3f(1., 1., lic_slice_z); glVertex2f((float)NMESH, (float)NMESH);
+    glTexCoord3f(0., 1., lic_slice_z); glVertex2f((float)0.,    (float)NMESH);
 
     /*
@@ -1222 +1332 @@
 
 
-void activate_one_volume_shader(int n_actual_slices){
+void activate_one_volume_shader(int volume_index, int n_actual_slices){
 
   cgGLSetStateMatrixParameter(m_mvp_vert_std_param, CG_GL_MODELVIEW_PROJECTION_MATRIX, CG_GL_MATRIX_IDENTITY);
@@ -1230 +1340 @@
 
   cgGLSetStateMatrixParameter(m_mvi_one_volume_param, CG_GL_MODELVIEW_MATRIX, CG_GL_MATRIX_INVERSE);
+  cgGLSetTextureParameter(m_vol_one_volume_param, volume[volume_index]->id);
   cgGLEnableTextureParameter(m_vol_one_volume_param);
-  cgGLSetParameter4f(m_render_param_one_volume_param, n_actual_slices, 0., 0., 0.);
+
+  //hack!
+  //if parameter.y == 0 : volume :0
+  //if parameter.y == 1 : volume :1
+  if(volume_index==0)
+    cgGLSetParameter4f(m_render_param_one_volume_param, n_actual_slices, 0., 0., 0.);
+  else if(volume_index==1)
+    cgGLSetParameter4f(m_render_param_one_volume_param, n_actual_slices, 1., 0., 0.);
   cgGLBindProgram(m_one_volume_fprog);
   cgGLEnableProfile(CG_PROFILE_FP30);
@@ -1241 +1359 @@
   cgGLDisableProfile(CG_PROFILE_FP30);
 
-  cgGLDisableTextureParameter(m_vel_tex_param);
-  cgGLDisableTextureParameter(m_pos_tex_param);
+  cgGLEnableTextureParameter(m_vol_one_volume_param);
+  //cgGLDisableTextureParameter(m_vel_tex_param);
+  //cgGLDisableTextureParameter(m_pos_tex_param);
 }
 
@@ -1364 +1483 @@
     */
     
-    activate_one_volume_shader(n_actual_slices);
+    activate_one_volume_shader(volume_index, n_actual_slices);
     glPopMatrix();
 
@@ -1500 +1619 @@
 #endif
 
+void draw_axis(){
+
+  glDisable(GL_TEXTURE_2D);
+  glEnable(GL_DEPTH_TEST);
+
+  //red x
+  glColor3f(1,0,0);
+  glBegin(GL_LINES);
+    glVertex3f(0,0,0);
+    glVertex3f(1.5,0,0);
+  glEnd();
+
+
+  //blue y
+  glColor3f(0,0,1);
+  glBegin(GL_LINES);
+    glVertex3f(0,0,0);
+    glVertex3f(0,1.5,0);
+  glEnd();
+
+
+  //green z
+  glColor3f(0,1,0);
+  glBegin(GL_LINES);
+    glVertex3f(0,0,0);
+    glVertex3f(0,0,1.5);
+  glEnd();
+
+  glEnable(GL_TEXTURE_2D);
+  glDisable(GL_DEPTH_TEST);
+}
+
+
 
 /*----------------------------------------------------*/
@@ -1511 +1663 @@
 
    //convolve
-   //lic();
+   lic();
 
    /*
@@ -1527 +1679 @@
    
    //advect particles
-   //psys->advect();
+   psys->advect();
 
    final_fbo_capture();
@@ -1547 +1699 @@
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
-   //glTranslatef(live_obj_x, live_obj_y, live_obj_z);
-   glTranslatef(0, 0, live_obj_z);
+   glTranslatef(live_obj_x, live_obj_y, live_obj_z);
 
    glRotated(live_rot_x, 1., 0., 0.);
@@ -1554 +1705 @@
    glRotated(live_rot_z, 0., 0., 1.);
 
-   /*
+   
+   glPushMatrix();
+
+   glScaled(volume[0]->aspect_ratio_width, 
+          volume[0]->aspect_ratio_height, 
+          volume[0]->aspect_ratio_depth);
+
+   
    //draw line integral convolution quad
    glBegin(GL_QUADS);
-   glTexCoord2f(0, 0); glVertex3f(0, 0, slice_z);
-   glTexCoord2f(1, 0); glVertex3f(1, 0, slice_z);
-   glTexCoord2f(1, 1); glVertex3f(1, 1, slice_z);
-   glTexCoord2f(0, 1); glVertex3f(0, 1, slice_z);
+   glTexCoord2f(0, 0); glVertex3f(0, 0, lic_slice_z);
+   glTexCoord2f(1, 0); glVertex3f(1, 0, lic_slice_z);
+   glTexCoord2f(1, 1); glVertex3f(1, 1, lic_slice_z);
+   glTexCoord2f(0, 1); glVertex3f(0, 1, lic_slice_z);
    glEnd();
-   */
+   
+   glPopMatrix();
+   
 
    //soft_display_verts();
-   //psys->display_vertices();
-
-   //render multiple volumes
+   
+   glPushMatrix();
+
+   glScaled(volume[0]->aspect_ratio_width,
+          volume[0]->aspect_ratio_height, 
+          volume[0]->aspect_ratio_depth);
+
+   psys->display_vertices();
+
+   glPopMatrix();
+
+   //render volume :0 
    //volume[0]->location =Vector3(0.,0.,0.);
    //render_volume(0, 256);
 
-   //render another but shifted using the same texture
-   volume[0]->location =Vector3(-0.5,-0.5,-0.5);
+   //render volume :1
+   volume[1]->location =Vector3(0., 0., 0.);
+   render_volume(1, 256);
 
    perf->enable();
-   render_volume(0, 256);
    perf->disable();
    fprintf(stderr, "pixels: %d\n", perf->get_pixel_count());
    perf->reset();
+
+   draw_axis();
 
    glDisable(GL_DEPTH_TEST);
@@ -1652 +1823 @@
                 break;
         case '+':
-                slice_z+=0.1;
+                lic_slice_z+=0.05;
                 get_slice_vectors();
                 break;
         case '-':
-                slice_z-=0.1;
+                lic_slice_z-=0.05;
                 get_slice_vectors();
                 break;
@@ -1667 +1838 @@
         case '3':
                 psys_x-=0.1;
+                break;
+        case 'w': //zoom out
+                live_obj_z-=0.1;
+                break;
+        case 's': //zoom in
+                live_obj_z+=0.1;
+                break;
+        case 'a': //left
+                live_obj_x-=0.1;
+                break;
+        case 'd': //right
+                live_obj_x+=0.1;
+                break;
+        case 'i':
+                init_particles();
                 break;
         default:

trunk/gui/vizservers/nanovis/nanovis.h

@@ -87 +87 @@
 float live_rot_z = 0.;
 
-float live_obj_x = -0.5;        //object translation location from the origin
-float live_obj_y = -0.5;
+float live_obj_x = -0.0;        //object translation location from the origin
+float live_obj_y = -0.0;
 float live_obj_z = -2.5;
 
@@ -95 +95 @@
 bool right_down = false;
 
-float slice_x=0, slice_y=0, slice_z=0.3;        //image based flow visualization slice location
+float lic_slice_x=0, lic_slice_y=0, lic_slice_z=0.3;//image based flow visualization slice location
 
 int win_width = NPIX;                   //size of the render window

trunk/gui/vizservers/nanovis/shaders/one_volume.cg

@@ -25 +25 @@
   float4 tex_coord = mul(modelViewInv, IN.TexCoord);
         
-#if 0
-  //1 component
-  float sample = tex3D(volume, tex_coord.xyz);
+  if(renderParameters.y==1){
+    //1 component
+    float sample = tex3D(volume, tex_coord.xyz);
 
-  if (sample.x < 0.) {
+    if (sample.x < 0.) {
       OUT.Color = float4(0.0, 0.0, 0.0, 0.0);
-  } else {
-      OUT.Color = float4(sample, sample, sample, sample/renderParameters.x);
+    } 
+    else{
+      //shades of red //more opaque
+      OUT.Color = float4(sample, 0, 0, 7*sample/renderParameters.x);
+
+      if(sample<0.6 && sample>0.4){
+        OUT.Color = float4(0, 0, 1, 10*sample/renderParameters.x);
+      }
+    }
   }
-#endif
 
+  else if(renderParameters.y==0){
 
-#if 1
-  //3 component
-  float3 sample = tex3D(volume, tex_coord.xyz).xyz;
-  OUT.Color = float4(3*sample, 2*length(sample)/renderParameters.x);
-#endif 
+    //3 component
+    float3 sample = tex3D(volume, tex_coord.xyz).xyz;
+    
+    OUT.Color = float4(3*sample, 2*length(sample)/renderParameters.x);
+  }
 
 //debug

trunk/gui/vizservers/nanovis/shaders/render_vel.cg

@@ -27 +27 @@
     ret.y = plane_normal.y * vel.y;
     ret.z = plane_normal.z * vel.z;
-
-    return float4(ret, 1);
+    
+    float3 tmp = ret*2.0 - float3(1,1,1);
+    return float4(tmp, 1);
 }

trunk/gui/vizservers/nanovis/shaders/update_pos.cg

@@ -18 +18 @@
             uniform samplerRECT pos_tex,
             uniform sampler3D vel_tex,
-            uniform float timestep
+            uniform float timestep,
+            uniform float3 scale
             ) : COLOR
 {
@@ -28 +29 @@
 
     float4 vel = float4(tex3D(vel_tex, pos.xyz).xyz, 0);
+
+    //reconstruct negative value
+    vel = (vel - float4(0.5,0.5,0.5,0.))*2.;
+    vel.x *= scale.x;
+    vel.y *= scale.y;
+    vel.z *= scale.z;
     ret = pos + vel*timestep;
-
+    //ret = pos;
 
 /*
+    //not drawing if the particle is out of bound
     if(ret.x<0 || ret.x>1 || ret.y<0 || ret.y>1 || ret.z<0 || ret.z>1 || time==0)
       ret = float4(1,0,0,0); 
     else
-      ret.w = time-1;
+      ret.w = 1;
 */
 
-    ret.w = 1.0;
+    ret.w = 1.;
     return ret;
 }