PDA

View Full Version : off topic - Volume Rendering question



Kal Mar
09-19-2011, 02:37 PM
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

mobeen
09-19-2011, 10:55 PM
Hi Kal,

the points are plotted in diffuse form.
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)
http://www3.ntu.edu.sg/home2007/mova0002/test/2DTF.png

Kal Mar
09-20-2011, 06:09 AM
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

mobeen
09-20-2011, 07:38 AM
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.

Kal Mar
09-20-2011, 12:05 PM
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

mobeen
09-20-2011, 08:06 PM
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.

Kal Mar
09-21-2011, 10:37 AM
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

mobeen
09-21-2011, 10:43 PM
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.

Kal Mar
09-22-2011, 05:10 AM
I'm sorry, my error was unintentional. Mobeen.

Sincerely,

Kal

Kal Mar
09-22-2011, 05:12 AM
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( &amp;x, &amp;y, &amp;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.

mobeen
09-22-2011, 04:41 PM
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])); }
}

Kal Mar
09-23-2011, 07:04 AM
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

mobeen
09-23-2011, 07:26 AM
So does it work?

Kal Mar
09-23-2011, 08:53 AM
Unfortunately, no.

mobeen
09-23-2011, 09:20 AM
Hmm then u need to check your code. I have already shared my code which produced the picture on the first page.

Kal Mar
09-23-2011, 02:13 PM
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

mobeen
09-24-2011, 06:31 AM
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.

Kal Mar
09-25-2011, 02:07 PM
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

mobeen
09-25-2011, 08:04 PM
I am glad it worked out fine for u.

Regards,
Mobeen

Kal Mar
10-05-2011, 05:07 AM
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

snowbin
11-15-2011, 07:34 PM
HI&amp;#65292; Kal Mar ,mobeen
I need your help, I am doing the work just as your guys did. my code is almost the same as yours. I just couldn't get the result like your guys show .
I stored the log-scale histogram into a 2D texture and then display that texture.the result is show in the follow picture.

glTexImage2D(GL_TEXTURE_2D, 0 , GL_RGBA8, 256, 221, 0, GL_RGBA, GL_UNSIGNED_BYTE, histogram);

glBegin(GL_QUADS);
glTexCoord2f(1, 1); glVertex2f(0, 0);
glTexCoord2f(0, 1); glVertex2f(1.0, 0);
glTexCoord2f(0, 0); glVertex2f(1.0, 1.0);
glTexCoord2f(1, 0); glVertex2f(0, 1.0);
glEnd();

thinks
Bin

mobeen
11-15-2011, 11:20 PM
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
Hi Kalmar,
Just use a larger histogram to accomodate the 16-bit values. and then setup the other params as required.

mobeen
11-15-2011, 11:20 PM
HI&amp;#65292; Kal Mar ,mobeen
I need your help, I am doing the work just as your guys did. my code is almost the same as yours. I just couldn't get the result like your guys show .
I stored the log-scale histogram into a 2D texture and then display that texture.the result is show in the follow picture.

glTexImage2D(GL_TEXTURE_2D, 0 , GL_RGBA8, 256, 221, 0, GL_RGBA, GL_UNSIGNED_BYTE, histogram);

glBegin(GL_QUADS);
glTexCoord2f(1, 1); glVertex2f(0, 0);
glTexCoord2f(0, 1); glVertex2f(1.0, 0);
glTexCoord2f(0, 0); glVertex2f(1.0, 1.0);
glTexCoord2f(1, 0); glVertex2f(0, 1.0);
glEnd();

thinks
Bin
Check your email.

snowbin
11-18-2011, 01:37 AM
Thanks, Mobeen.
I got the packet you send me.

mobeen
11-18-2011, 07:15 AM
Hmm so can u now see the gradient magnitde histogram?

snowbin
12-04-2011, 01:06 AM
yeah! I wrote my own code already.

Thanks a lot