off topic - Volume Rendering question

Hi guys

Please, I need help, I will be grateful. Any help is welcome.

I am working in volume rendering and I have a trouble with the 2D Transfer Function (especifically displaying the scalar data and their gradient magnitude).

I computed the gradient (central difference) and find the gradient magnitude, but when I show the gradient magnitude as points with OPENGL over the transfer function domain (as in the paper of Kniss02 - Mult. TF), I don´t visualize the archs between materials as I should want. Instead, the points are plotted in diffuse form.

I revised my alg, I change to Log-Scale, but it is useless.

Thanks,

Kal

Hi Kal,

Could u show us what u r getting right now. Usually the plot should be something like this (this is for visible male head dataset)

Hi Mobbem

You are right. I don’t understand how to restrict the gradient value for get the arch. Please give me some ideas

My plot of ct foot dataset is showed in:

http://www.megaupload.com/?d=594IJKV4

(I’m sorry, I could not add a image to this reply)

Thanks,

Kal

Hi Kal,
No this is not what u should get. It seems to me that the grad mag. histogram is incorrect. How are u calculating it at the moment?
Giving u some ideas. Like the normal image histogram, the log scaled gradient magnitude is used instead. Since ur dataset (volume) is 3d, you would move over each slice, calculate the grad mag., log scale it and then increment the histogram bin corresponding to the gradient magn. Each slice will add to the same 2D histogram.

Try it out and let me know if u cannot do it. I will share my implementation but u need to show me some progress.

Hi Mobeen

I calculated the grad mag(its are right, I finish check it) by central difference,

Vector3 sample1, sample2; //typedef vector (x,y,z)

for (int z = 0; z < dim[Z]; z++) {
for (int y = 0; y < dim[Y]; y++) {
for (int x = 0; x < dim[X]; x++) {
sample1.x = SampleVol8(data, x - 1, y, z);
sample2.x = SampleVol8(data, x + 1, y, z);
sample1.y = SampleVol8(data, x, y - 1, z);
sample2.y = SampleVol8(data, x, y + 1, z);
sample1.z = SampleVol8(data, x, y, z - 1);
sample2.z = SampleVol8(data, x, y, z + 1);
vec = vec_sub(sample2, sample1); //difference
vec = vec_mul(0.5, vec); //div by 2.0
gradMag[index] = vec_mag(vec); //magnitud
index++;
}
}
}

but I, at the moment, don’t perform any compute of histogram (simply, i just plot scalar versus grad mag), the result was the plot that I send.

for (int z = 0; z < dim[Z]; z++) {
for (int y = 0; y < dim[Y]; y++) {
for (int x = 0; x < dim[X]; x++) {
//histogram graphics
posx = data[x + (y * dim[X]) + (z *dim[X]
* dim[Y])];
posy = gradMag[index];
glBegin(GL_POINTS);
glColor3f(0.5, 0.5, 0.5);
glVertex3f( posx, posy, 0.0f);
glEnd();
index++;
}
}
}

Now, I am thinking how to “increment” it, as you said. What will be the axis X, will be the scalar data ? and the axis y will be frequency of grad mag ? i’m sorry, don’t understand it. By instance, a scalar 20 in foot data set has 2836 ocurrences and has very different grad mag, ex. 0.477121 (3), 0.539591 (2), 0.615224 (7), 0.80103(5)…, so on.

You mean in some way accumulate it and plot this values?

Thanks,

Kal

Hi Kal,
Well it is a histogram of gradient magnitudes which contains the intensities (the voxel value) on X axis and gradient magnitudes on the Y axis. For a simple image histogram u do this



//Assuming 8-bit greyscale image with intensities from 0-255
//H -> image height
//W -> image width

int histo[256]; //histogram
for y:0 to H
   for x:0 to W
      intensity = GetPixel(x,y)
      histo[intensity]++;
   end for
end for

In our case since, we have a 3D volume, we need to do this,


//Assuming 8-bit dataset
//XDIM = x dimension of dataset
//YDIM = y dimension of dataset
//ZDIM = z dimension of dataset
//pVolume is the pointer to volume data
//our histogram is 2D now with voxel intensity on X-axis and gradient magnitude on Y-axis.
int histo[256][256];
for z:0 to ZDIM
   for y:0 to YDIM
      for x:0 to XDIM
         index = ((z*YDIM)+y)*XDIM+x;
         voxel = pVolume[index];
         grad  = GetGradient(x,y,z);
         gradMag = sqrt(grad.x*grad.x+ grad.y*grad.y + grad.z*grad.z);
         histo[gradMag][voxel]++;  
      end for
   end for 
end for

//log scale histogram
for j=0 to 256
   for i=0 to 256
      histo[j][i] = log(histo[j][i]);
   end
end for

Do note that

  1. the gradient is not normalized in the above code.
  2. the log value will be returned in floating point so u would need to scale/cast it based on the output type u have.

Another thing I would like to point out. For rendering this, it is better to use a 2D texture rather than GL_POINTS. Just pass the histogram to glTexImage2D (u have to take care of flipping the Y values and rescale the data based on the texture size). If for testing to render as points, you should put glBegin(GL_POINTS) and glEnd() outside the nested loops.


glBegin(GL_POINTS);
for (...)
   for (...)
      for (...)
         glVertex2f(...);
      end for
   end for
end for
glEnd();

See if this helps.

Hi Mobeem

Your help served me and was very useful.

But excuse me, I’ve trouble yet; after log-scale hist2D (without scaling), I plot the points and it has a bell form, similar to hist1D.

Maybe the rounding are changing the precision (I am using integer por mapping coordinates), but the fact is that I can not still see the arches.

My plot was:

Thank,

Kal

Hi Kal,
My name is M O B E E N not MOBEEM.
Second for you case, could u show us how u r creating your histogram? Leave the log scaling for the moment just tell us how r u doing it currently.

I’m sorry, my error was unintentional. Mobeen.

Sincerely,

Kal

Hi Mobeen

I post mi code for compute the histogram:

  
//8-bit dataset

//dim[X] dimension of dataset in X
//dim[Y] dimension of dataset in Y
//dim[Z] dimension of dataset in Z

uint32 gradMag,scalar;
vector3 grad; //vector(x,y,z)
 
for (int z = 0; z < dim[Z]; z++) {
 for (int y = 0; y < dim[Y]; y++) {
  for (int x = 0; x < dim[X]; x++) {
	scalar = SampleVolume(x, y, z);
	grad = GetGradient(x, y, z);
	gradMag = (uint32) vector_magnitude(grad);
	histogram2D[gradMag][scalar]++; //compute 2D histogram		
  } 
 } 
}

histogram2D was intialized with:

 
uint32 histogram2D[256][256];

And one sample from volume is getting with:


uint8 SampleVolume(int x, int y, int z) {
   ClampIndexVolume( &x, &y, &z);
   return (uint8) dataset[x + (y * dim[X]) + (z * dim[X]* dim[Y])];
} 

The method GetGradient is:


vector3 GetGradient( int x, int y, int z) {
	vector3 vecGrad, sample1, sample2;
	sample1.x = SampleVolume( x - 1, y, z);
	sample2.x = SampleVolume( x + 1, y, z);
	sample1.y = SampleVolume( x, y - 1, z);
	sample2.y = SampleVolume( x, y + 1, z);
	sample1.z = SampleVolume( x, y, z - 1);
	sample2.z = SampleVolume( x, y, z + 1);
	vecGrad = vector_sub(sample2, sample1);
	vecGrad = vector_div(0.5, vecGrad); 
	return vecGrad;
} 

when I add the number of samples in the histogram2D, effectively corresponds to the number of data elements.

Hi Kal,
I think this line

vecGrad = vector_div(0.5, vecGrad); 

should be

vecGrad = vector_div(2, vecGrad); 

. Is it possible for u to share ur dataset so that I can tell u what output u should get may be the dataset contains the arches that u r getting right now?
I am doing exactly the same thing as u r doing here. This is how i do it.


int ih[256][256]={0};
for(int k=0;k<ZDIM;k++)
{
   int zInd1 = (k==0)?k:k-1;
   int zInd2 = (k==ZDIM-1)?k:k+1;
    
   for(int j=0;j<YDIM;j++)
   {
      int yInd1 = (j==0)?j:j-1;
      int yInd2 = (j==YDIM-1)?j:j+1;
      for(int i=0;i<XDIM;i++)
      {				
         int xInd1 = (i==0)?0:(i-1);
	 int xInd2 = (i==XDIM-1)?i:i+1;

	 float gradient[3];
	 int vox = ((k*YDIM + j)*XDIM + i);
	 int val  = pVolume[vox];
	 gradient[0]=(pVolume[(k*YDIM     + j)    *XDIM +xInd2]-pVolume[(    k*YDIM+j)    *XDIM+xInd1])*0.5f;
	 gradient[1]=(pVolume[(k*YDIM     + yInd2)*XDIM +i    ]-pVolume[(    k*YDIM+yInd1)*XDIM+i    ])*0.5f;
	 gradient[2]=(pVolume[(zInd2*YDIM + j)    *XDIM +i    ]-pVolume[(zInd1*YDIM+j)    *XDIM+i    ])*0.5f;				
	 int gradMag =  (sqrt(gradient[0]*gradient[0]+gradient[1]*gradient[1]+gradient[2]*gradient[2]));
 
	 ih[gradMag][val]++;
      }
   } 
}
//Log scale the histogram
for(int j=0;j<256;j++) {
   for(int i=0; i< 256; ++i){
	ih[j][i] = (float) log(float(ih[j][i]));    }
}

Thanks for your answer Mobeen,

I corrected the error in the division. Very thanks.

Currently I am testing with 8-bist dataset downloaded from volvis.org, especifically the foot and engine dataset. I know that these have arches. The maximum and minimun grad-mag (without root square) for the foot and engine I find, (Min= 0 Max= 154),(Min= 0 Max= 129), respectively. With that I could discard if the compute of gradien is correct.

About your work, what are you doing?, I begin my work of volume segmentation, but as you see, I have many troubles.

Have a nice day!

Kal

So does it work?

Unfortunately, no.

Hmm then u need to check your code. I have already shared my code which produced the picture on the first page.

Hi Mobeen,

I checked that my code is right computing the gradient histogram, but I don’t get the expected result.

Please, how are you mapping the scatterplots?

My mapping is:


glBegin(GL_POINTS);
for (int grad_mag = 0; grad_mag < 256; grad_mag++) {
  for (int scalar = 0; scalar < 256; scalar++) {
	glVertex2f(scalar, hist2D[grad_mag][scalar]);
  }
}
glEnd();

Thanks,

Kal

OK. For rendering of histogram, I first normalize the histogram by dividing with the max and min values. Then scale the y values based on the size of the histogram on screen.
So roughly it is this way.


int min=1000, max=-1000;
//find the min and max value of histogram
for(int j=0;j<256;j++) {
   for(int i=0;i<256;i++) {
      if(histo[j][i]>max)
        max = histo[j][i];
      if(histo[j][i]<min)
        min = histo[j][i];
   }
}

//float normHisto[256][256]; // in 0-1 range
//now scale the histogram
for(int j=0;j<256;j++) {
   for(int i=0;i<256;i++) {
      normHisto[j][i] = ((float)histo[j][i])/(max-min);
   }
}

//scale the histogram for display
for(int j=0;j<256;j++) {
   for(int i=0;i<256;i++) {
 //scale value and store in an array
         }
}

I suggest u store the scaled histogram into a 2D texture and then display that texture.

Hi Mobeen,

I am very grateful for your valuable help. The problem was my code that was not correctly reading the data.

I send the images that I get from the head and foot data set

Have a good day.

Regards,

Kal

I am glad it worked out fine for u.

Regards,
Mobeen

Hi Mobeen

Excuse me, maybe do you have any suggestions to render the hist2D as 2D texture for 16-bit datasets ?

Thank you in advance.

Have a nice day,

Kal