ray casting problem

Dear all I have a ray casting problem, I must perform raycasting in a local sphere S at a voxel x , but I don't know how to increment d or t in this formule, p = t d + p0 can any one give me a complete code please? thanks in advance

[Bere]

I mean how to produce the first and next rays?

please help me.

Answer 1

[Grzegorz Knor]

Are you familiar with this file?

http://www.mathworks.de/matlabcentral/fileexchange/26852-a-fast-voxel-traversal-algorithm-for-ray-tracing

[Grzegorz Knor]

BTW:

http://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

Answer 2

[Bere]

it was so useful, but I want to implement a code that perform raycasting in a local sphere S at a voxel x with radius R. A total number of NS rays are uniformly distributed cast from the center x to the sphere’s surface. For each ray, the casting stops when it reaches the boundary of the sphere or when it encounters the boundary of a local structure. Can you help me how to change this code to do this? .

[Walter Roberson]

Your first step is to figure out how to do a uniform random distribution over a sphere. It is easy to get that step wrong. I suggest you search the Newsgroup for back postings by Roger Stafford as he has posted the algorithm a few times.

What do you want to *do* with the rays when they stop? If your answer is "draw them", then be advised that you are going to have difficulty representing a uniform random 3-space distribution to a fundamentally pixelized display.

What is your representation of the local structures?

Answer 3

[Bere]

I want to measure boundary by comparing the image intensity at the current position p on the ray and at the spherical center x. When I(p)-I(x)/I(x)>e , this indicates a local boundary (e is the threshold that controls the sensitivity of local boundary detection). After casting in all directions, only the rays whose lengths are longer than 0.8 R are preserved. please help me. thanks alot

Answer 4

[Bere]

please help me ,I have a project and I must implement it.

[Walter Roberson]

Which formula are you using for the uniform random distribution over the sphere?

What is your representation of the local structures whose presence must be detected?

Answer 5

[Bere]

I don't have any special formula, just a uniform distribution over a ray from a center of sphere, it must detect in coronary vessels a tube, a tube junction, a plane and a sphere. In the tube and the tube junction, the rays are oriented towards the tube directions. While, in the plane and the sphere,the rays are scattered without giving any explicit direction. thanks alot

[Walter Roberson]

If you do not have a formula for the uniform distribution over the sphere, then you have *not* done the research I suggested in to Roger Stafford's postings on the subject.

When you do not do your part in working on *your* question, it leaves people uninterested in helping you.

Read this:

http://www.mathworks.com/matlabcentral/newsreader/view_thread/150991

Perhaps someone else will clarify the difference between what an object _is_ and how you _represent_ the object.

Answer 6

[Bere]

thanks alot for your reply,

I use this code for uniform distribution,

I test all of them and think it's good for 3d images,

function [x,y,z] = rsphere(N)
z = 2 * rand(1,N) - 1;
rho = sqrt(1-z.^2);
phi = pi * (2 * rand(1,N) - 1);
x = rho .* cos(phi);
y = rho .* sin(phi);

can anybody help me how to trace a ray between two points in center and surface of the sphere ?

thanks in advance

[Walter Roberson]

Now you need to tell us about the representation of the various objects you have. What data structure do you use to hold the information about them.

For example, the tube: you must have some kind of mathematical model of it, or some kind of approximation of the surface of it. If I were to specify a particular coordinate (x1,y1,z1) then how would it be possible to tell whether that coordinate was inside the tube, at the surface of the tube, or outside the tube?

If a point were determined to be at the surface of the tube, then how would it be possible to tell what the normal is at that point, so that the angle of reflection can be determined?

Is the tube made out of glass? With uniform thickness? With uniform coloration? Is it transparent or translucent or completely opaque? If it is transparent or translucent then how do you represent the index of refraction so that it would be possible to model the bending of light by the tube?

The tube junction: with the information you have given so far, we might guess that the tube is made out of glass or plastic, but if so then it is very unlikely to be a constant thickness all through the junction; and especially if it is glass, chances are strong that it has different indices of refraction at different points on the surface and as you go in to the object. How is the information about that represented?

Answer 7

[Bere]

When I(p)-I(x)/I(x)>e , this indicates a local boundary ,I just want to trace the intensity of voxels in a ray and save those rays which have above condition, The distribution of these rays describes the geometric features of the local shape, which is then measured to discriminate between the false step-edge responses of cardiac cavities and the desired responses of coronary arteries. To quantify the distribution of the rays, connected component labeling on the sphere is performed to cluster the rays oriented in the same local direction into M groups. thank you.

[Bere]

I extracts the vessels before and they are gray scale with low resolution. it must discriminate

between the false step-edge responses of cardiac

cavities and the desired responses of coronary

arteries.

[Walter Roberson]

What are I(p) and I(x), and e ? And is your formula perhaps incorrect? You have I(x)/I(x) which would simplify to 1 when I(x) is non-zero.

For a given (x0,y0,z0) point, how do we determine I(p) and I(x) ? As you do not mention y or z, is the implication that I is independent of y and z?

Are you perhaps trying to indicate that what you have is a 3D array of intensities and you want to do the ray tracing with respect to that intensity array? That you do not have a mathematical model of the tube and so on, just an array of intensities?

[Bere]

ofcourse I have pixel dimension, it is [0.36 , 0.36,0.4] .

Answer 8

[Bere]

I(p) is intensity of centeral voxel and I(x) is intensity of current voxel in ray,

yes I must paranthesize the formulea , it is (I(p)-I(x))/I(x)>e e is a treshhold.

 no, x is a indices of voxels, I have 3d images in 512*512*272 dimension, and I don't know how I must specify a voxel in it.

yes I have just an array of intensity and then I must specify the type of structure by another formulea.

thanks again.

[Walter Roberson]

Okay, I think I finally understand. I will have to think a bit about a suitable algorithm.

Should the voxels be treated as cubes each of uniform intensity within the cube, with a hard boundary on each? Or should interpolation be done?

Which voxel is your central voxel ? As your dimensions are all even, the center of the array is at the point of contact of 8 voxels. And to cross-check: the rays are to have their origin at the center of the centeral voxel, right?

[Bere]

yes, a linear interpolation must be done,

The raycasting parameters R and e varied with the optimal scale

F(x) and R(x) = 3(F(x) + 1), F(x) varies between 1-3.

sorry the vosel dimensions are:[0.32 0.32 0.4] and size of image after downsampling must be 256*256*272.

would you please tell me how you compute 8 pixels around the central voxel?

central voxel varies in different parts of image, yes the rays have their origin at the center of the central voxel.

thank you.

[Walter Roberson]

Ah, I did not notice before that you are doing local ray casting from a number of different points, and I had not realized before the significance of the surface of the sphere.

I take it that R(x) is the radius (shouldn't that be R(p) if p is your current central voxel?) Does F(x) have any significance to the ray casting, or is it just "something that is needed to calculate the radius" as far as the ray casting goes?

Having the voxels not be equal length in each dimension adds some challenge, but should be workable.

This might be easier (faster) if no linear interpolation is done.

Oh my... potentially up to 27 voxels involved in a single interpolation. This I'm going to have to think on more carefully.

Answer 9

[Bere]

yes, it is R(p), and F(x) has no significance to the ray casting, it is a parameter in the filter in previous step and has no relation to ray casting. thank you.

[Bere]

would you please tell me how you calculate the number of voxels?

[Walter Roberson]

Consider a point right at the center of a 3x3x3 cube. You might perhaps define the value for that point to be the value associated with the center voxel, a single value.

Now move the point. As it approaches the outside of the center voxel, your need for linear interpolation requires that the values of adjacent voxels be taken in to account in proportion to the distance of the point to the center of that voxel.

Suppose the point were to be moved to an exact vertex of the central voxel. Should only the voxels along the 3 axes (up/down, left/right, front/back) together with the center be taken in to account for the interpolation? Logically, No: at the vertex, the distance to the center of the center voxel is the same as the distance to the center of voxel that the diagonal continuation of the direction of the center towards the point, so if the contribution relative to distance is what is taken in to account, the diagonal needs to be weighted the same as the center.

Likewise as you move around near the edges of the central voxel at other locations (e.g., near one of the edges), you clearly need to take in to account the voxels in some of the other directions including diagonals on other planes.

As the cube is 3x3x3, that is up to 27 different voxels that might be involved in the calculation.

What is the maximum number of voxels that need to be considered for any one of the calculations? Well, you might have *defined* the value at the very center of the central voxel to be the same as the value of the center voxel, but as soon as you "wiggle" the point position even a little off of the center, that definition cannot hold for a linear interpolation, so the the number is going to depend on exactly how you define your linear interpolation. If you cull the voxels according to which half-planes the point under consideration is in, because the point might have only wiggled on one of the three axes, you might only cull a single 3x3 face, leaving 18 voxels to take in to account. Even then it is hard to justify culling a voxel that is distance 1/2 + epsilon away while still justifying using the contribution of the central voxel at distances further than 1/2 away from the center (i.e., beyond the outside of the maximum sphere that could be inscribed in the cube). Culling by half-plane becomes even more suspect when you take in to account that your voxels are not actually cubes.

Linear interpolation is based upon contribution of value proportional to distances, and on a simple straight line does ignore the value on "the other side", but I need to think more on how to mathematically formulate that cutoff in a manner that is fair for cuboids. Or, I suppose, I could just use interpn() and hope that it takes care of the details fairly.

Perhaps linear interpolation should not be used if the point falls within the inscribed sphere? Note that the inscribed sphere would have the radius equal to the smallest voxel dimension, not touching the walls along the larger dimension.

Answer 10

[Fareena]

Dear Bere

Did you solve this? If yes; please provide me with the code.