Google Earth -> LDraw/Minecraft ("Voxel" Decomp
» >
» Goto page Previous  1, 2, 3, 4, 5, 6, 7, 8, 9, 10  Next
» View previous topic :: View next topic  
oh okay.

What are you taking a determinant of that is yielding zero? That sounds like a degenerate triangle and/or some inappropriately parallel lines. Are you sure the model contains only tris and no quads?
<- Please Support on Lego Ideas
Yes, all quads. Here's the code in question. tri is a 3x3 array, rows are points, columns are [x, y, z]. Given a known x and y or x and z or y and z, it finds the missing point on the triangle, if possible.


Code:
    def geom_cart2bary(self,tri,x,y,z):
        P = array([x,y,z])
        if x == None:
            a = 1
            b = 2
            c = 0
        elif y == None:
            a = 0
            b = 2
            c = 1
        else:
            a = 0
            b = 1
            c = 2

        det = (tri[1,b]-tri[2,b])*(tri[0,a]-tri[2,a]) + \
              (tri[2,a]-tri[1,a])*(tri[0,b]-tri[2,b])

        if det == 0:        # Incoming div-by-0
            print("Determinant: Div-by-0 warning")
            return [0., 0., 0., x, y, z]

        # Compute barycentric coordinates
        l1 = ((tri[1,b]-tri[2,b])*(P[a]-tri[2,a]) + \
              (tri[2,a]-tri[1,a])*(P[b]-tri[2,b]))/det;
        l2 = ((tri[2,b]-tri[0,b])*(P[a]-tri[2,a]) + \
              (tri[0,a]-tri[2,a])*(P[b]-tri[2,b]))/det;
        l3 = 1. - l1 - l2

        rval = [l1, l2, l3, x, y, z]
        rval[3+c] = (l1*tri[0,c]) + (l2*tri[1,c]) + (l3*tri[2,c])

        return rval
Geopipe
can you give some example inputs producing success and failure modes?

Just to clarify: tri is a matrix of row vectors defining a triangle, and [x,y,z] are a barycentric coordinate triplet missing a coordinate?
<- Please Support on Lego Ideas
Tri is a matrix of row vectors defining a triangle. [x, y, z] are cartesian coordinates missing one of the coordinates. The desired output is [l1, l2, l3, x, y, z], where x, y, and z are all filled, and l1, l2, and l3 are the barycentric lambda coordinates.

Sample success matrix:
Code:
[[  84.68746953   59.83120068  155.70715938]
 [  82.88576245   60.64850894  152.60306357]
 [  83.41658276   58.7063019   152.60306357]]


Sample failure matrix:
Code:
[[ 49.42879067  88.10337266  56.30088594]
 [ 57.97150024  83.39647856  56.89353193]
 [ 49.42879067  88.10337266  56.89353193]]


Notice that [0,0] = [2,0] and [0,1] = [2,1] (and because the z-axis span is the smallest, z will be the missing coordinate). This will make (tri[1,b]-tri[2,b])*(tri[0,a]-tri[2,a]) + (tri[2,a]-tri[1,a])*(tri[0,b]-tri[2,b]) of course be zero. And that will be true even if the triangle points get shuffled. What to do.

Edit: Unless I shuffle rows 1 and 2 in this case, I believe... That would give:
Code:
[[ 49.42879067  88.10337266  56.30088594]
 [ 49.42879067  88.10337266  56.89353193]
 [ 57.97150024  83.39647856  56.89353193]
]
For which we get the det = (nonzero)*(nonzero) + (nonzero)*(nonzero). And I guess I can then just reshuffle the lambdas at the end.

Edit #2: Nope, although that makes the terms nonzero, it also makes the sum cancel out to zero. I should have expected it wasn't that easy.
Geopipe
What matrix are you trying to take the determinant of? Because that is almost certainly not the formula for a 3x3 determinant.

WabbitEmu agrees with me: I just plugged those in time to my calculator to do some determinant calculations and got nonzero results for the original, or the rowswapped version. Then I looked more closely at the formula and it doesn't appear to have the correct terms unless I'm really bad at mental algebra (it foils to 8 terms, not 6).

In this case, |tri| = -583.9326763.

[edit]

As a side note, the matrices would absolutely be singular if you swap only a[0,2] and a[2,2]. Are all of the triangles where you're having trouble right triangles where the 90° bend is along a primary axis?

[edit 2]

What happens if you replace your code with this:


Code:
    def geom_cart2bary(self,tri,x,y,z):
        P = array([x,y,z])
        if x == None:
            a = 1
            b = 2
            c = 0
        elif y == None:
            a = 0
            b = 2
            c = 1
        else:
            a = 0
            b = 1
            c = 2

        det = numpy.lingalg.det(tri[:,array([a,b,c])])

        if det == 0:        # Incoming div-by-0
            print("Determinant: Div-by-0 warning")
            return [0., 0., 0., x, y, z]

        # Compute barycentric coordinates
        l1 = ((tri[1,b]-tri[2,b])*(P[a]-tri[2,a]) + \
              (tri[2,a]-tri[1,a])*(P[b]-tri[2,b]))/det;
        l2 = ((tri[2,b]-tri[0,b])*(P[a]-tri[2,a]) + \
              (tri[0,a]-tri[2,a])*(P[b]-tri[2,b]))/det;
        l3 = 1. - l1 - l2

        rval = [l1, l2, l3, x, y, z]
        rval[3+c] = (l1*tri[0,c]) + (l2*tri[1,c]) + (l3*tri[2,c])

        return rval

Also, how do you choose which columns to swap? I don't understand the progression from 1,2,0 to 0,2,1 to 0,1,2 ( I would have expected that middle list to be perhaps 2,0,1). Also also, swapping columns only effects the sign of the determinant since it's an alternating multilinear function (by rows or columns), so that makes me think you're calculating a determinant of some kind even if it's not actually the determinant of |tri| (it was scale-invariant under row swapping).

[edit 2]
Oh, I see, I think you're trying to take the determinant of


But I'm pretty sure you can't do that, because the barycentric coordinate formulation is given in terms of x and y only, so you need to calculate planar coordinates instead of just dropping a coordinate.

Let's adjust your matrix in the following manner:


Code:

# fiddle with this a bit more to choose the axises of rotation
# and proper indexing if/when
# the missing coordinate/smallest edge is not the third one
P = matrix([[x],[y],[z]])


trimatrix = matrix(tris).transpose()
offset = trimatrix[:,0]
centered = tris - offset # translate to origin
theta_x = atan2(centered[2,2],centered[1,2])
rx = matrix([[1,0,0],[0,cos(-theta_x),-sin(-theta_x)],[0,sin(-theta_x),cos(-theta_x)]])
intermediate1 = rx*centered
theta_z = atan2(intermediate1[1,2],intermediate1[0,2])
rz = matrix([[cos(-theta_z),-sin(-theta_z),0],[sin(-theta_z),cos(-theta_z),0],[0,0,1]])
theta_x2 = atan2(intermediate2[2,1],intermediate2[1,1])
rx2 = matrix([[1,0,0],[0,cos(-theta_x2),-sin(-theta_x2)],[0,sin(-theta_x2),cos(-theta_x2)]])
planar = rx2*intermediate2

Pplanar = rx2*rz*rx*(P-offset)

# Do something here to set tolerances
# and take care of rounding error, if it seems like it matters

T = matrix([[(planar[0,0]-planar[0,2]),(planar[0,1]-planar[0,2])],[(planar[1,0]-planar[1,2]),(planar[1,1]-planar[1,2])]])
det = numpy.linalg.det(T)
if det == 0:        # Incoming div-by-0
    print("Determinant: Div-by-0 warning")
    return [0., 0., 0., x, y, z]
lam12 = inv(T) * (Pplanar[:2]-planar[:2,2])
lambda1 = lam12[0]
lambda2 = lam12[1]
lambda3 = 1 - lambda1 - lambda2
rval = [lambda1,lambda2,lambda3,x,y,lambda1*trimatrix[2,0]+lambda2*trimatrix[2,1]+lambda3*trimatrix[2,2]]
<- Please Support on Lego Ideas
I'm computing the determinant of the T matrix as per this Wikipedia article. The problems do indeed seem to arise when there's a 90-degree bend along a primary axis. And I tried swapping rows, not columns.
Geopipe
I think if you mix the last codeblock I posted above with your values of a,b,c, and my reindexing trick it should work.


Code:
        if x == None:
            a = 1
            b = 2
            c = 0
        elif y == None:
            a = 0
            b = 2
            c = 1
        else:
            a = 0
            b = 1
            c = 2
reindex = array([a,b,c])
trimatrix = matrix(tris[reindex]).transpose()
#instead of trimatrix = matrix(tris).transpose()
P = matrix([[x],[y],[z]])[reindex]
#instead of P =matrix([[x],[y],[z]])
<- Please Support on Lego Ideas
Will this take into account the fact that one of the variables is unknown? i.e., will it mind that one of x, y, and z is None, which in the P = line will turn one of those elements into a NaN?
Geopipe
that's what the reindexing in the previous codesnippet takes care of. It always puts the missing one in last place.
<- Please Support on Lego Ideas
Gotcha, that makes sense. So the first problem I encounter with this (other than the missing intermediate2 that we worked out) is at Pplanar = rx2*rz*rx*(P-offset), where it complains about subtracting a float from a None. I resolved that by just setting P[2] to 0.0 all the time (in the hopes that since it's "unknown", its value doesn't matter), but I'm relatively sure that's not the correct solution. After that "fix", I get so many triangles the converter claims come out empty that I'm pretty sure something is wrong, but I also eventually get "TypeError: unsupported operand type(s) for -: 'NoneType' and 'numpy.int64'" in an operation that indicates lambda1*trimatrix[2,0]+lambda2*trimatrix[2,1]+lambda3*trimatrix[2,2] is coming out as None.

Also, doesn't the rval line require a remapping to reverse the coordinate remapping at the beginning?
Geopipe
hmmm, I'm afraid we need a proxy value to carry through in place of Z, and then solve some kind of equation at the end.

Let's do an experiment. Create this class:

Code:

class Unknown(object):
   def __init__(self,name,sterm=0,mterm=1):
      self.name = name
      self.reciprocal = False
      self.sterm=sterm
      self.mterm=mterm
      
   def __repr__(self):
      return "%s((%s*Unknown(%s))+sterm)" % ("1/" if self.reciprocal else "",repr(self.mterm),repr(self.name),repr(self.sterm))
   
   def __str__(self):
      return "((%s*%s)+sterm)" % ("1/" if self.reciprocal else "",str(self.mterm),str(self.name),str(self.sterm))
   
   def __mul__(self, other):
      self.sterm = self.sterm * other
      self.mterm = self.mterm * other
      return self
   
   def __div__(self, other):
      self.sterm = self.sterm / other
      self.mterm = self.mterm / other
      return self
   
   def __add__(self, other):
      self.sterm = self.sterm + other
      return self
   
   def __sub__(self, other):
      self.sterm = self.sterm - other
      return self
      
   def __rmul__(self, other):
      return self.__mul__(other)
   
   # Division is not commutative!
   def __rdiv__(self, other):
      self.reciprocal = not self.reciprocal
      return self.__div__(1/other)
   
   def __radd__(self, other):
      return self.__add__(other)
   
   # Subtraction is not commutative!
   def __rsub__(self, other):
      return self.__sub__(-other)
   
   def __neg__(self, other):
      return self.__mul__(-1)

Then do this after the initial setup

P[2] = Unknown("Z")

and then let fly, and print out the lambdas at the end and see what it gives you. I recommend using test data, not real data/not trying to write data output Wink

w.r.t the coordinate remapping, the lambdas should be spitting out in exactly the same fashion you had before, the x/y/z probably do need your
Code:
rval[3+c]
trick instead of how I set it up..
<- Please Support on Lego Ideas
okay, so I've produced a patched Unknown class, and as far as I can tell, it works just fine with my numpy. I also noticed that the line which says
Code:
centered = tris - offset
should be
trimatrix instead of tris


Code:
class Unknown(object):
   def __init__(self,name,sterm=0,mterm=1,reciprocal=False):
      self.name = name
      self.reciprocal = reciprocal
      self.sterm=sterm
      self.mterm=mterm
      
   def __repr__(self):
      return "%s((%s*Unknown(%s))+%s)" % ("1./" if self.reciprocal else "",repr(self.mterm),repr(self.name),repr(self.sterm))
   
   def __str__(self):
      return "%s((%s*%s)+%s)" % ("1./" if self.reciprocal else "",str(self.mterm),str(self.name),str(self.sterm))
   
   def __mul__(self, other):
      oterm = (1./other) if self.reciprocal else other
      return Unknown(name=self.name,sterm=self.sterm*oterm,mterm=self.mterm*oterm,reciprocal=self.reciprocal)
   
   def __div__(self, other):
      return self.__mul__(1./other)
   
   def __add__(self, other):
      if self.reciprocal:
         # This is where it gets hairy
         denom = Unknown(name=self.name,sterm=self.sterm,mterm=self.mterm,reciprocal=True)
         return Unknown(name=self.name,sterm=(self.sterm*other+1)*denom,mterm=(self.mterm*other)*denom,reciprocal=False)
      else:
         return Unknown(name=self.name,sterm=self.sterm+other,mterm=self.mterm,reciprocal=False)
      self.sterm = self.sterm + other
      return self
   
   def __sub__(self, other):
      return self.__add__(-other)
      
   def __rmul__(self, other):
      return self.__mul__(other)
   
   # Division is not commutative!
   def __rdiv__(self, other):
      return Unknown(name=self.name,sterm=self.sterm,mterm=self.mterm,reciprocal = not self.reciprocal).__rmul__(other)
   
   def __radd__(self, other):
      return self.__add__(other)
   
   def __rsub__(self, other):
      return self.__neg__().__radd__(other)
   
   def __neg__(self):
      return Unknown(name=self.name,sterm=-self.sterm,mterm=-self.mterm,reciprocal=self.reciprocal)
   
   def __call__(self, value):
      term = self.sterm + (self.mterm * value)
      return 1/term if self.reciprocal else term


The only tricky part is getting it to solve for itself, and that shouldn't actually be too hard since everything is linear.
<- Please Support on Lego Ideas
Here's the result of my unit test:

Code:
>>> t = kmz2mc.Tri2Voxel(None)
>>> t.geom_prep([0, 0, 0],[100, 100, 100])
Reserving 101 x 101 x 101 array...
>>> tri
array([[ 93.7878877 ,  30.0793981 ,  56.30088594],
       [ 76.68651914,  34.28190747,  56.30088594],
       [ 78.89287129,  38.28630181,  56.30088594]])
>>> t.geom_cart2bary(tri, 80., 35., None)
[((-0.0*Unknown('Z'))+((-0.0*Unknown('Z'))+0.15027257070615613)), ((-0.0*Unknown('Z'))+((0.0*Unknown('Z'))+0.5126940849587779)), ((0.0*Unknown('Z'))+((0.0*Unknown('Z'))+((0.0*Unknown('Z'))+((-0.0*Unknown('Z'))+0.33703334433506604)))), 80.0, 35.0, ((-0*Unknown('Z'))+((-0*Unknown('Z'))+((-0*Unknown('Z'))+((0.0*Unknown('Z'))+((0.0*Unknown('Z'))+((0.0*Unknown('Z'))+((0.0*Unknown('Z'))+((-0*Unknown('Z'))+56.300885940000001))))))))]
Clearly the last two lines shouldn't be:
Code:
        rval = [lambda1,lambda2,lambda3,x,y,z]
        rval[3+c] = lambda1*trimatrix[2,0]+lambda2*trimatrix[2,1]+lambda3*trimatrix[2,2]
Geopipe
Looks good (last 2 lines included)! Just need to solve for dem Zs! I'll try to make a quick solver Smile
<- Please Support on Lego Ideas
elfprince13 wrote:
Looks good (last 2 lines included)! Just need to solve for dem Zs! I'll try to make a quick solver Smile
Yeah, I'm currently pondering it myself; just wanted to point out that everything is looking good. Smile I'm also trying to wrap my brain around this custom object with custom operators idea; this is new to me.

Edit: looks like lambda1-3 can at least just be lambda1.sterm.sterm..., but I guess I need an operator or method (float()?) that "simplifies" these Unknown objects, discarding the 0*name mterms.
Geopipe
I have a simplifier set up. Do you have scipy as well as numpy?
<- Please Support on Lego Ideas
elfprince13 wrote:
I have a simplifier set up. Do you have scipy as well as numpy?
I do indeed; I even have it loaded, because pymclevel requires it.
Geopipe
Meet me back here in 15 minutes Smile
<- Please Support on Lego Ideas
elfprince13 wrote:
Meet me back here in 15 minutes Smile
Will do! I'll investigate the thin-triangle issue in the meantime.
Geopipe
The solver works pretty well, though it will probably work better if you can pass in a reasonable guess at z.

Haven't tested the last couple lines with the returns - should work, but if it complains about broadcasting the function call with **try_solve in the arguments, you may have to make those assignments one at a time.



Code:
import scipy.optimize as scop
class Unknown(object):
   def __init__(self,name,sterm=0,mterm=1,reciprocal=False):
      self.name = name
      self.reciprocal = reciprocal
      self.sterm=sterm
      self.mterm=mterm
      
   def __repr__(self):
      return "%s((%s*Unknown(%s))+%s)" % ("1./" if self.reciprocal else "",repr(self.mterm),repr(self.name),repr(self.sterm))
   
   def __str__(self):
      return "%s((%s*%s)+%s)" % ("1./" if self.reciprocal else "",str(self.mterm),str(self.name),str(self.sterm))
   
   def __mul__(self, other):
      oterm = (1./other) if self.reciprocal else other
      return Unknown(name=self.name,sterm=self.sterm*oterm,mterm=self.mterm*oterm,reciprocal=self.reciprocal)
   
   def __div__(self, other):
      return self.__mul__(1./other)
   
   def __add__(self, other):
      if self.reciprocal:
         # This is where it gets hairy
         denom = Unknown(name=self.name,sterm=self.sterm,mterm=self.mterm,reciprocal=True)
         return Unknown(name=self.name,sterm=(self.sterm*other+1)*denom,mterm=(self.mterm*other)*denom,reciprocal=False)
      else:
         return Unknown(name=self.name,sterm=self.sterm+other,mterm=self.mterm,reciprocal=False)
      self.sterm = self.sterm + other
      return self
   
   def __sub__(self, other):
      return self.__add__(-other)
      
   def __rmul__(self, other):
      return self.__mul__(other)
   
   # Division is not commutative!
   def __rdiv__(self, other):
      return Unknown(name=self.name,sterm=self.sterm,mterm=self.mterm,reciprocal = not self.reciprocal).__rmul__(other)
   
   def __radd__(self, other):
      return self.__add__(other)
   
   def __rsub__(self, other):
      return self.__neg__().__radd__(other)
   
   def __neg__(self):
      return Unknown(name=self.name,sterm=-self.sterm,mterm=-self.mterm,reciprocal=self.reciprocal)
   
   #Evaluate one level deep
   def __call__(self, **variables):
      tterm = Unknown(self.name)
      mdict = {}
      sdict = {}
      for n,v in variables.iteritems():
         if self.name == n:
            tterm = v
         if type(self.mterm) == Unknown and n in self.mterm:
            mdict[n] = v
         if type(self.sterm) == Unknown and n in self.sterm:
            sdict[n] = v
      term = (self.sterm(**sdict) if sdict else self.sterm) + (tterm * (self.mterm(**mdict) if mdict else self.mterm))
      return 1/term if self.reciprocal else term
   
   def __contains__(self, key):
      ret = False
      if self.name == key:
         ret = True
      if type(self.mterm) == Unknown:
         ret |= self.mterm.__contains__(key)
      if type(self.sterm) == Unknown:
         ret |= self.sterm.__contains__(key)
      return ret
   
   def findaroot(self, **guesses):
      args = sorted(guesses.keys())
      init_vals = array([guesses[a] for a in args])
      def redict(argorder,values):
         return {a:v for a,v in zip(args,values)}
      
      def solvable(array_values):
         return self.__call__(**redict(args,array_values))
         
      def sstatus(x,f):
         print "Trying %s to get %f" % (str(x),self.__call__(**redict(args,x)))
      
         
      result = scop.root(solvable,init_vals,method='excitingmixing',tol=1e-10)
      if result.success:
         ret = redict(args,result.x)
      else:
         ret = {}
      return ret





if x == None:
   a = 1
   b = 2
   c = 0
elif y == None:
   a = 0
   b = 2
   c = 1
else:
   a = 0
   b = 1
   c = 2 

reindex = array([a,b,c])
trimatrix = matrix(tris[reindex]).transpose() 
P = matrix([[x],[y],[z]])[reindex]
P[2]=Unknown('z')

offset = trimatrix[:,0]
centered = trimatrix - offset # translate to origin
theta_x = atan2(centered[2,2],centered[1,2])
rx = matrix([[1,0,0],[0,cos(-theta_x),-sin(-theta_x)],[0,sin(-theta_x),cos(-theta_x)]])
intermediate1 = rx*centered
theta_z = atan2(intermediate1[1,2],intermediate1[0,2])
rz = matrix([[cos(-theta_z),-sin(-theta_z),0],[sin(-theta_z),cos(-theta_z),0],[0,0,1]])
intermediate2=rz*intermediate1
theta_x2 = atan2(intermediate2[2,1],intermediate2[1,1])
rx2 = matrix([[1,0,0],[0,cos(-theta_x2),-sin(-theta_x2)],[0,sin(-theta_x2),cos(-theta_x2)]])
planar = rx2*intermediate2

Pplanar = rx2*rz*rx*(P-offset)

# Do something here to set tolerances
# and take care of rounding error, if it seems like it matters

T = matrix([[(planar[0,0]-planar[0,2]),(planar[0,1]-planar[0,2])],[(planar[1,0]-planar[1,2]),(planar[1,1]-planar[1,2])]])
det = numpy.linalg.det(T)
if det == 0:        # Incoming div-by-0
    print("Determinant: Div-by-0 warning")
    return [0., 0., 0., x, y, z]

lam12 = inv(T) * (Pplanar[:2]-planar[:2,2])
lambda1 = lam12[0,0]
lambda2 = lam12[1,0]
lambda3 = 1 - lambda1 - lambda2
rval = [lambda1,lambda2,lambda3,x,y,lambda1*trimatrix[2,0]+lambda2*trimatrix[2,1]+lambda3*trimatrix[2,2]]

rval = array([lambda1,lambda2,lambda3,x,y,z])
rval[3+c] = lambda1*trimatrix[2,0]+lambda2*trimatrix[2,1]+lambda3*trimatrix[2,2]

# feel free to put better guesses if you have them
# but it's linear so I think it should be okay?
try_solve = (rval[3+c] - P[2,0]).findaroot(z=0)
if try_solve:
    unsolved_indices = array([1,2,3,3+c])
    rval[unsolved_indices] = rval[unsolved_indices](**try_solve)
else:
    print "Found no solution to equation %s=0" % (str(rval[3+c] - P[2]))
    return [0., 0., 0., x, y, z]

return rval
<- Please Support on Lego Ideas
  
Register to Join the Conversation
Have your own thoughts to add to this or any other topic? Want to ask a question, offer a suggestion, share your own programs and projects, upload a file to the file archives, get help with calculator and computer programming, or simply chat with like-minded coders and tech and calculator enthusiasts via the site-wide AJAX SAX widget? Registration for a free Cemetech account only takes a minute.

» Go to Registration page
Google Earth -> LDraw/Minecraft ("Voxel" Decomp
» Goto page Previous  1, 2, 3, 4, 5, 6, 7, 8, 9, 10  Next
» View previous topic :: View next topic  
Page 2 of 10
» All times are UTC - 5 Hours
 
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum

 

Advertisement