Query geometry in specific coordinate of 3D space

[mrwolfy]

If I have a piece of geomety in a scene, how much can I find out about the geometry at a particular area in coordinate space?

What I want to do specifically is to devide my scene into 3D cubic grid units, the area inside a cube defined by 4 points in xyz space. I want to query weather or not geometry exhists in each cube. You can see in my illustration I have drawn a single arbitrary "unit" of grid space as I am envisioning it. I want to find out if geometry from my object intesects this "cube area".

I also would like to be able to know the color of the object in that location. Ideally it would be an average color, the average color of the texture that intersects the area of the 3D grid unit.

Is this possibel in Maya?

Any insight would be very much appreciated!

[nowayfra]

I guess this goes into that direction PIXELS by PATRICK JEAN does. Practically what you could do with standard nodes is:

1. use polyCubes for your 3D cubic grid. You then always know where the inside and outside volume ends, depending on surface normals. You will need this later.
2. You will need to create a nearestPointOnMesh node (i now call this thing "npom") for your mesh objects in your scene and a closestPointOnSurface for all nurbsSurfaces (both of them you get with "createNode nearestPointOnMesh/closestPointOnSurface")
3. To query for example if a mesh object intersects one of your poly grid cubes, feed the npom two times:

First: 
sceneMeshShape.outMesh -> npom.inMesh
gridCube.translate -> npom.inPosition
(value of grid cube translate should be identical with grid cube pivot, but at least must be located inside the gridCube) If you are a hardcore mathguy you now could do some angular stuff with the vector from grid cubes center to the nearest point of your scene mesh and your grid cube edges and their vectors to find out wether the objects intersect or not, but if not: save the result position (the nearest point on mesh vector value, npom.position) as resultPoint1. Now reconnect the npom:

Second: 
gridCubeShape.outMesh -> npom.inMesh
resultPoint1 -> npom.inPosition
Now finally interesting are those two vectors: the vector(1) FROM nearest point on scene mesh TO nearest point on gridCube* AND the vector(2) of the nearest point on your gridCube NORMAL (npom.normal). If your nearest scene mesh point lies inside the gridCube volume: the ANGLE between vector(1) and vector(2) is less than 90 degrees (because the vector targets from inside the gridCube to outside like the cridCube normals do. If your nearest point on scene mesh lies outside the grid cube, the angle is greater than 90 degrees.

This is how to work with mesh objects. Nurbs surfaces need the closestPointOnSurface but the workflow will be the same. Now the corresponding color value. Simply use the resulting paramaterU and parameterV values of your npom node and feed those values into a "colorAtPoint" command, its source color is the assigned materials input color node (file, ramp or something). You can travel the connections from your secene mesh to your color node with the "listConnections" command. Analyze the Haypershade connections for this. If you use some more complex shader networks for your meshes color value, bake them first. The colorAtPoint doesnt query cumulated colors at a single color node. To get an average color you could shoot some random "rays" from your gridCube pivot into the direction of your closest point of scene mesh and qery if they are inside, then get the color again.
You of course will need some basic scripting knowledge for this, but you surely didnt expect to find a button for this.
Good luck.

*this vector value is the result of gridCubes nearest point MINUS secene mesh nearest point

[joojaa]

everybody (the user, not every one) has written this script its called legoize. But he hasn't made it public but teh entire workflow is described nonetheless in a older post.  Anyway one could ask for the script instead of reimplementing it (yeah i know where one can get it but its not mine so i can't distribute it). I dont think it does color but anyway.

Thing is the nearest point of mesh is a bit inefficent method. theres a way faster method, which really comes into play whan you make a hundredthousand to million cubes, on acomplex mesh.

Getting avreage color is a bit hard tough, since theres really no way to derivate over the area. So you must prefilter the data somehow.

Anyway one could take a existing quadtree maker and modify it for this need also if ones really industrious

[mrwolfy]

Thanks guys!

nowayfra, thanks for your detailed instructions, wow! I do have some scripting and am familiar with colorAtPoint, but it will take me a while to absorb all this good information nonetheless.

JooJaa, No problem about the average color, I think a sampled point of color anywhere in the "area"/"volume" is fine. Any more specifics would be fantastic help. What is "a way faster method"?

[mrwolfy]

Hi Nowayfra,

I think I have everything working up to the point where I am supposed to measure the angle between vectors. I tried using the angleBetween command, which returns four numbers, three with the euler flag. I am not sure what to do with these values. I am studying on it, but any insight would be helpful.

Thanks![url]../../../users/nowayfra/forums[/url]

[nowayfra]

Hi,
the angleBetween is for getting an angle and axis you need to rotate one vector to match the direction of another vector. Hence those many return values. To simply get the angle between two vectors, use "angle" command, for example  angle vector1 vector2. It then returns the angle in radians so you could convert them with rad_to_deg for better readability but you just can use the radiant value. Then an angle less than 90° =  is less than 1.570796 (pi / 2), greater than 90° = greater than 1.570796.

[joojaa]

Or even shorter is to take the dot product between the vectors and if its bigger than 0 they point towards the same locality. This saves you one call of arcsin and 2 calls of normalisation. Anyway.

If you actually intend to do many cubes Insetad of just one here and one there then read this before you continue:

The real problem with this method is that its a O(n^3) method, what everybody did was order the problem a bit better to get down to a manageable O(n^2). Just to point out what this means. It means that if you take a

1x1x1 = 1 operation
10x10x10 = 1000 operations
100x100x100 = 1000000 operations
1000x 1000 x 1000 = 1000000000 operations
10000x10000x10000 = 1000000000000 operations

Now how long does one opearton take, lets take a very optimistic view on the speed of a computation to 10 milliseconds or 0.01 seconds but as you can see this doesnt really matter. This means the 10 cubed takes 10 seconds 100 cubed takes 10000 seconds or 166 minutes or 2.7 hours where the 1000 cubed is starting to be a bit tedious 2777.777 hours or ~116 days, now the 10000 cubed takes 317 YEARS! So its safe to say that this method will never be used over 100x100x100 cubes no matter how slow you can push the individual computation. Speed doenst really mater offcourse even if you did a 100 times faster youd only manage 1000 cubed at max at very uncofortable speeds.

If on the otherhand the computation was  O(n^2) then youd be down a notch making 1000x1000x1000 as speedy as 100x100x100 and making 100 cubed ~17 seconds versus 2.7 hours. However all this said this CAN be made into a very manageable O(n*log(n)) operation quite easily or a longwindedO(n) operation perhaps if ones very smart a stochastic O(log(n)).

THIS IS WHAT PROGRAMMING SKILL IS ABOUT, versus mere scripting skill.

Just so wa stand on the same page. Sometimes it pans out to think how much data you are processing.

[nowayfra]

Yes of course we surely are aware what cubic growth means. And yes bringing this up to api level including a quadtree algorithm for volume exclusion is a good idea (besides this is also a good idea to speed mel scripts up). But i am really curious about a different and "way faster" approach to query if volumes intersect.

[joojaa]

Well thing is ratracing allready does this. So for example each and every time you start mentalray with raytracing, then mentalray does just this. It takes each polygon and splits it up into a smallish cube, with which it can do this. And it does so quite fast in order of seconds per million of polygons.

Since you know in advance the size of each cube you can just, split up each triangle into each cube it passes trough. If yoru really good you can do this in a single level of operations. Hence youd get a o(n) But seriously you probably get something like O(n*log(n)) if you account for testing what cubes you will hit. But usually thats not what O notation is about.

Theres a really easy way to cheat. Shoot rays from one direction, and record every intersection along that ray and then populate a array (if you hit at voxel 3 and voxel 5 and voxerl 12 and voxel 102  you know voxels 3-5 and 12-102 are solid). Once your array is populated find the outlying volume of the cubes by only including the voxels that dont have a neigbour thats included. Now since your shooting only form one side you get a O(n^2) solution because your not processing one side (also you can select the smallest side to do the raycasting).

Method 2 is to actually let maya dice your geometry with poly split, and then bound each diced part. This also is faster, but destroys the mesh in the process. Again its bringing you down in O range but not much.

Also theres a built in fucntion as of maya 2008 or something where you can use the fill function of nDynamics to do the same thing at full api speed. And done by a proper programmer too (Duncan), so i presume its quite efficient

Now you can expand this to actually find the volume but in practice this is enough for most meshes. You can shoot rays in mel if you compile the rayitersect plugin in devkit.

[mrwolfy]

@Nowayfra

Hi, I've set this up, per your instructions. Obviously I did something wrong. I'm still looking into it, but if it is clearly obvious to you what I have missed, would you mind pointing it out to me?

Thanks!!!

//create nearestPointOnMesh nodes
createNode nearestPointOnMesh -name npom1;
createNode nearestPointOnMesh -name npom2;
string $theCubes[] = `ls -type transform "cube\_*"`;
for($obj in $theCubes)
{
//make connections as described in the instructions by nowayfra: http://tinyurl.com/2av2323
//part 1, First
connectAttr -f sourceGeometryShape.outMesh npom1.inMesh;
connectAttr -f ($obj + ".translate") npom1.inPosition;
//store the results of part 1, First
float $vector1f[] = getAttr("npom1.position");
vector $vector1 = getAttr("npom1.position");
//part 2, Second
connectAttr -f ($obj + ".outMesh") npom2.inMesh;
setAttr ($obj + ".resultPoint1X") $vector1f[0];
setAttr ($obj + ".resultPoint1Y") $vector1f[1];
setAttr ($obj + ".resultPoint1Z") $vector1f[2];
connectAttr ($obj + ".resultPoint1X") npom2.inPositionX;
connectAttr ($obj + ".resultPoint1Y") npom2.inPositionY;
connectAttr ($obj + ".resultPoint1Z") npom2.inPositionZ;
//store the results of part 2, Second
vector $vector2 = getAttr("npom2.normal");
float $angle = `angle $vector1 $vector2`;
//disconnect everything
disconnectAttr  sourceGeometryShape.outMesh npom1.inMesh;
disconnectAttr ($obj + ".translate") npom1.inPosition;
disconnectAttr  ($obj + ".outMesh") npom2.inMesh;
disconnectAttr ($obj + ".resultPoint1X") npom2.inPositionX;
disconnectAttr ($obj + ".resultPoint1Y") npom2.inPositionY;
disconnectAttr ($obj + ".resultPoint1Z") npom2.inPositionZ;
if ( $angle > 0.5) {
    delete $obj;
}
}

fyi, I am storing the value of resultPoint1 in a vector attribute that I added to each cube.

Before

After

The python to create the cube setup.

import maya.cmds as mc
\# create a block grid of cubes
mult =  3;
count = 0
cubes = []
y = mult - mult
while y <= mult *2:
    x = mult * -1
    while x <= mult:
        z = mult * -1
        while z <= mult:
            cube = mc.polyCube( name = ('cube\_' + str(count) ) )
            ## add a vector attribute to each cube to store closestPointOnMesh result.
            mc.addAttr( cube[0], ln = 'resultPoint1', at = 'double3' )
            mc.addAttr( cube[0], ln = 'resultPoint1X', at = 'double', p = 'resultPoint1' )
            mc.addAttr( cube[0], ln = 'resultPoint1Y', at = 'double', p = 'resultPoint1' )
            mc.addAttr( cube[0], ln = 'resultPoint1Z', at = 'double', p = 'resultPoint1' )
            mc.setAttr( cube[0] + '.resultPoint1', 0,0,0, type = 'double3')
            mc.setAttr( cube[0] + '.resultPoint1', keyable = True)
            mc.setAttr( cube[0] + '.resultPoint1X', keyable = True)
            mc.setAttr( cube[0] + '.resultPoint1Y', keyable = True)
            mc.setAttr( cube[0] + '.resultPoint1Z', keyable = True)
            count +=1
            cubeName=cube[0]
            cubes.append(cubeName)
            mc.move(x, y, z)
            z += 1
        x +=1
    y +=1
\##create source geometry
mc.polySphere(name = 'sourceGeometry')
mc.move(0,3,0)
mc.setAttr('sourceGeometry.scaleX', 3)
mc.setAttr('sourceGeometry.scaleY', 3)
mc.setAttr('sourceGeometry.scaleZ', 3)

[nowayfra]

Ok, i have to admit that the method described is a bit inaccurate... the nearest point on scene mesh may lie outside the gridCube but the cube can be intersected nonetheless. You could fix this with some tricky additional lines, but i guess this would end in an approximation only. Anyway, try this:
First freeze transformations on your scene mesh, then also on your gridCubes. The nearestPointOnMesh seem to return values in object space, so interaction in world space doesnt work properly. Also you didnt use the right angle value for comparison (the "angle" command gives its value in radiants, means 90° equals value pi/2).
Just create your sphere and select its transform, then try those lines:

proc makeCubeGrid(int $cubeGridDiv)
{
    string $cube[];
    string $cubes[] = {};
    string $sel[] = ls -long -sl;
    float $bBox[] = polyEvaluate -b -ae $sel[0];

    float $w = ($bBox[1] - $bBox[0]) / $cubeGridDiv;
    float $h = ($bBox[3] - $bBox[2]) / $cubeGridDiv;
    float $d = ($bBox[5] - $bBox[4]) / $cubeGridDiv;

    float $p[3];
    for($i=0;$i<$cubeGridDiv;$i++)
    {
        $p[0] = $bBox[0] + (($bBox[1] - $bBox[0]) / $cubeGridDiv) * $i;
        for($k=0;$k<$cubeGridDiv;$k++)
        {
            $p[1] = $bBox[2] + (($bBox[3] - $bBox[2]) / $cubeGridDiv) * $k;
            for($n=0;$n<$cubeGridDiv;$n++)
            {
                $p[2] = $bBox[4] + (($bBox[5] - $bBox[4]) / $cubeGridDiv) * $n;
                $cube = polyCube -w $w -h $h -d $d -sx 1 -sy 1 -sz 1 -ax 0 0 0 -ch 0;
                move ($p[0] + $w / 2) ($p[1] + $h / 2) ($p[2] + $d / 2) $cube[0];
                makeIdentity -apply true -t 1 -r 1 -s 1 -n 0 $cube[0];
                $cubes[size $cubes] = $cube[0];
            };
        };
    };    
    string $group = group -em -n "cubes";
    parent $cubes $group;
};

makeCubeGrid(10);

//now select the speheres transform again and execute this:

string $sel[] = ls -long -sl;
makeIdentity -apply true -t 1 -r 1 -s 1 -n 0 $sel[0];
string $cubes[] = listRelatives -children "cubes";
string $cubeShapes[] = listRelatives -allDescendents -type mesh "cubes";
string $npom1 = createNode nearestPointOnMesh;
string $npom2 = createNode nearestPointOnMesh;
$sel = listRelatives -children $sel[0];
connectAttr ($sel[0] + ".outMesh") ($npom1 + ".inMesh");
for($i=0;$i<size $cubes;$i++)
{
    float $t[3] = xform -q -ws -piv $cubes[$i];
    setAttr ($npom1 + ".inPosition") $t[0] $t[1] $t[2];
    vector $p1 = getAttr ($npom1 + ".position");
    connectAttr -f ($cubeShapes[$i] + ".outMesh") ($npom2 + ".inMesh");
    setAttr ($npom2 + ".inPosition") ($p1.x) ($p1.y) ($p1.z);
    vector $p2 = getAttr ($npom2 + ".position");
    vector $n = getAttr ($npom2 + ".normal");
    float $angle = angle(($p2 - $p1), $n);
    if($angle > 1.5707) setAttr ($cubes[$i] + ".visibility") 0;
};

if you want to measure time have a look at the "timerX" command

[joojaa]

wouldnt a naiveapproach of testing wether a vertex exists within bounding box suffice as long as we are talking of approximations. Because atleast for the example it does work. And thast helofalot faster.

I mean its not like hes interested in the interior volume!

[nowayfra]

Yes but consider a simple 1x1 plane put in those cubes. You then end with four cubes - one at each vertex.

[mrwolfy]

That works!

[joojaa]

Yes but consider a simple 1x1 plane put in those cubes. You then end with four cubes - one at each vertex.

I am aware of that. however. since tesselating is a relatively cheap methodology compared to the testof finding closest point. Especially since we possibly recreate the octree over and over again.

See we dont actually know the border constraints of the project. Sometimes the engineer in us tends to go for the solution that allways works. But sometimes the simplest thing that MIGHT work is the right thing to do. If its a throwaway code project then it doesnet matter. Since our solution doesnt really scale well anyway then bringing in a solution thats simpler but still works with a few manual gloss overs and DOES scale after that may be the thing to do. So theres a huge reisk that we are over engineering this.

I wanted to test doing the problem as a bucket sort with splitting of polygons by all buckets and have python api code in place that works, its slow but its a n log n slope magnitude solution, uses a lot fo memory tough so it should be using considerably less memory if converted to c++. And im currently multitherading the solution. For most meshes im currently doing this voxelisation process in near real time. (100*100*100 cubes with a mesh of 10,000 points at 20 fps, and i could be doing considerably better)