LoginSignup
5
1
HoudiniAdvent Calendar 2023Day 7
@kit2cuz(Kaz Kitaguchi)

交線を求める

Last updated at Posted at 2023-12-06

過去に交点を求める記事を書きましたが、今回は交線を求める記事を書きたいと思います。

Houdiniでフェースとフェースの交線を求める場合は、
Intersection Analysis SOPを使用して、そのノードのOutput Intersection Segmentsチェックボックスを有効にすることです。

IntersectionAnalysis

今回は、SOPを使用せずにWrangleで交線を求めたいと思います。
そういえば？VEX関数にintersect関数があるのでそれを使用すれば求まり・・・ません:cold_sweat:
そのintersect関数は光線(Ray)と面との交差ですね。
Wrangleで面と面の交線を求めるとなると独自に数式を求めてそれをコードにする必要があります。

では、数学的に面と面の交線はどうやって求めれば良いでしょう？
交線の方向はすぐに思いつくことでしょう。

交線の方向L = 面1の法線と面2の法線の両方に直交するベクトル。
つまり、面1の法線と面2の法線の外積で求まりますね。
外積

後は、交線が通る任意の1点Pを求めれば
交線 = P + α*L
という式ができますね。

平面1:

a1X+b1Y+c1Z=d1

平面2:

a2X+b2Y+c2Z =d2

とした場合、それぞれの平面を通る任意の1点は、X,Y,Zどれかを0として決め打ちすることで、この連立方程式を解くことができます。
交線方向LのX成分が0でない場合(YZ平面上にない場合)、Xを0と決め打ちすることでY,Zの値が求まります。
同様にLのY成分が0でない場合、Yを0と決め打ちすることでX,Zの値が求まり、
LのZ成分が0でない場合、Zを0と決め打ちすることでX,Yの値が求まります。

これをWrangleで実装すると

PrimitiveWrangle
vector nul = 0;
int success = 0;
vector curPrimNormal = prim_normal(0,0,nul);
vector curPrimPoint = primattrib(0,"P",0,success);
vector targetPrimNormal = prim_normal(1,0,nul);
vector targetPrimPoint = primattrib(1,"P",0,success);
//平面1と平面2の交線方向
vector intersectionDir = normalize(cross(curPrimNormal,targetPrimNormal));

//平面1の式の係数
float a1 = curPrimNormal.x;
float b1 = curPrimNormal.y;
float c1 = curPrimNormal.z;
float d1 = a1*curPrimPoint.x + b1*curPrimPoint.y + c1*curPrimPoint.z;
//平面2の式の係数
float a2 = targetPrimNormal.x;
float b2 = targetPrimNormal.y;
float c2 = targetPrimNormal.z;
float d2 = a2*targetPrimPoint.x + b2*targetPrimPoint.y + c2*targetPrimPoint.z;

vector p = 0;
float detA = 0.0;
if(intersectionDir[0]!=0.0){
    detA = b1*c2 - c1*b2;
    p = set(
        0.0,
        ( c2*d1 - c1*d2 )/detA,
        -( b2*d1 - b1*d2 )/detA);
}
else if(intersectionDir[1]!=0.0){
    detA = c1*a2 - a1*c2;
    p = set(
        -( c2*d1 - c1*d2 )/detA,
        0.0,
        ( a2*d1 - a1*d2 )/detA);
}
else if(intersectionDir[2]!=0.0){
    detA = a1*b2 - b1*a2;
    p = set(
        ( b2*d1 - b1*d2 )/detA,
        -( a2*d1 - a1*d2 )/detA,
        0.0);
}
int newLine = addprim(0, "polyline");
int newPtNum1 = addpoint(0, p-100.0*intersectionDir);
addvertex(0, newLine, newPtNum1);
int newPtNum2 = addpoint(0, p+100.0*intersectionDir);
addvertex(0, newLine, newPtNum2);

無限交線

これで、交線を求めることができました。
ただ、これは2枚の平面同士の無限交線の求め方で、実際に欲しいのは複数枚の平面同士の有限交線だと思います。
コードは以下のような感じで書きます。

PrimitiveWrangle
int isInside(vector refPos;int listPoints[]){
    int result = 0;
    int success = 0;
    if(len(listPoints)<3) return result;
    
    vector refPnt,pnt1,pnt2;
    int listTargetPoints[]=listPoints;
    push(listTargetPoints,listTargetPoints[0]);

    float originalArea = 0.0;
    refPnt = pointattrib(1,"P",listTargetPoints[0],success);
    for(int i =0;i<len(listTargetPoints)-3;i++){
        pnt1 = pointattrib(1,"P",listTargetPoints[i+1],success);
        pnt2 = pointattrib(1,"P",listTargetPoints[i+2],success);
        originalArea += 0.5*length(cross(pnt1-refPnt,pnt2-refPnt));
    }
    
    float currentArea = 0.0;
    refPnt = refPos;
    for(int i =0;i<len(listTargetPoints)-1;i++){
        pnt1 = pointattrib(1,"P",listTargetPoints[i],success);
        pnt2 = pointattrib(1,"P",listTargetPoints[i+1],success);
        float tmpArea = 0.5*length(cross(pnt1-refPnt,pnt2-refPnt));
        currentArea += tmpArea;
    }
    if( abs(currentArea - originalArea)<1e-6) result = 1;//微小面積対策として、精度厳しくする
    
    return result;
}

void getIntersection(vector A;vector B;vector C;vector D;int result;vector out)
{
    result = 1;
    if(A==B || C==D){
        result = 0;
        return;
    }
    vector n1 = normalize(B-A);
    vector n2 = normalize(D-C);
    if(abs(1.0-abs(dot(n1,n2)))<1e-6){//平行だった場合
        result = 0;
        return;
   }
    vector AC = C - A;
    float d1 = ( dot(n1,AC) - dot(n2,AC)*dot(n1,n2) ) / ( 1.0 - dot(n1,n2)*dot(n1,n2) );
    float d2 = ( -dot(n2,AC) + dot(n1,AC)*dot(n1,n2) ) / ( 1.0 - dot(n1,n2)*dot(n1,n2) );

    if( abs(d1)<1e-6 ) d1 = 0.0f;
    if( abs(d2)<1e-6 ) d2 = 0.0f;
    if(d1<0.0f || d2<0.0f) {
        result = 0;
        return;
    }
    
    if( abs(d1-distance(A,B)) <1e-6  ) d1 = distance(A,B);
    if( abs(d2-distance(C,D)) <1e-6  ) d2 = distance(C,D);
    
    if(d1>distance(A,B) || d2>distance(C,D)){
        result = 0;
        return;
    }
    vector p1 = A + d1*n1;
    if(distance(B,p1)<1e-6){
        result = 0;
        return;
    }
    out = p1;
    return;
}

void createIntersection(vector startPnt;vector endPnt;int targetPrimNum)
{
    int listPoints[] = primpoints(1,targetPrimNum);
    int newPntCount = 0;
    vector listNewPtPos[];
    if(isInside(startPnt,listPoints)==1){
        push(listNewPtPos,startPnt);
        newPntCount+=1;
    }
    if(isInside(endPnt,listPoints)==1){
        push(listNewPtPos,endPnt);
        newPntCount+=1;
    }
    if(newPntCount<2){
        push(listPoints,listPoints[0]);
        for(int i=0;i<len(listPoints)-1;i++){
            int success = 0;
            vector startPos = pointattrib(1,"P",listPoints[i],success);
            vector endPos = pointattrib(1,"P",listPoints[i+1],success);
            vector out;
            getIntersection(startPnt,endPnt,startPos,endPos,success,out);
            if(success==1){
                if(out!=startPnt && out!=endPnt){
                    push(listNewPtPos,out);
                    newPntCount+=1;
                    if(newPntCount==2) break;
                }
            }
        }
    }
    if(len(listNewPtPos)==2){
        int newLine = addprim(0, "polyline");
        
        foreach(vector eachNewPtPos;listNewPtPos){
            int newPtNum = addpoint(0, eachNewPtPos);
            addvertex(0, newLine, newPtNum);
        }
            
        addprimattrib(0,"line",0);
        setprimattrib(0,"line",newLine,1,"set");
    }
}

void getPrimIntersection(int curPrimNum;int targetPrimNum;int check;vector result1;vector result2)
{
    vector nul = 0;
    int success = 0;
    vector outOrigin=0;
    vector outVector=0; 
    vector curPrimNormal = prim_normal(0,curPrimNum,nul);
    vector curPrimPoint = primattrib(0,"P",curPrimNum,success);
    vector targetPrimNormal = prim_normal(1,targetPrimNum,nul);
    vector targetPrimPoint = primattrib(1,"P",targetPrimNum,success);
    vector intersectionDir = normalize(cross(curPrimNormal,targetPrimNormal));

    float a1 = curPrimNormal.x;
    float b1 = curPrimNormal.y;
    float c1 = curPrimNormal.z;
    float d1 = a1*curPrimPoint.x + b1*curPrimPoint.y + c1*curPrimPoint.z;
    float a2 = targetPrimNormal.x;
    float b2 = targetPrimNormal.y;
    float c2 = targetPrimNormal.z;
    float d2 = a2*targetPrimPoint.x + b2*targetPrimPoint.y + c2*targetPrimPoint.z;

    check = 0;
    vector p = 0;
    float detA = 0.0;
    if(intersectionDir[0]!=0.0){
        detA = b1*c2 - c1*b2;
        p = set(
            0.0,
            ( c2*d1 - c1*d2 )/detA,
            -( b2*d1 - b1*d2 )/detA
        );
        check = 1;
    }
    else if(intersectionDir[1]!=0.0){
        detA = c1*a2 - a1*c2;
        p = set(
            -( c2*d1 - c1*d2 )/detA,
            0.0,
            ( a2*d1 - a1*d2 )/detA
        );
        check = 1;
    }
    else if(intersectionDir[2]!=0.0){
        detA = a1*b2 - b1*a2;
        p = set(
            ( b2*d1 - b1*d2 )/detA,
            -( a2*d1 - a1*d2 )/detA,
            0.0
        );
        check = 1;
    }
    else{
        check = 0;
    }
    p = dot( curPrimPoint - p , intersectionDir ) * intersectionDir + p;
    if(check==1){
        outOrigin = p;
        outVector = intersectionDir;
    } 
    else{
        outOrigin = 0;
        outVector = 0;
    }
    if(check==1){
        int listPoints[] = primpoints(0,curPrimNum);
        vector startPnt = outOrigin - 100*outVector;
        vector endPnt = outOrigin + 100*outVector;
        vector results[]=array(startPnt,endPnt);
        int newPntCount = 0;
        push(listPoints,listPoints[0]);
        for(int i=0;i<len(listPoints)-1;i++){
            vector startPos = pointattrib(0,"P",listPoints[i],success);
            vector endPos = pointattrib(0,"P",listPoints[i+1],success);
            vector out;
            getIntersection(startPnt,endPnt,startPos,endPos,success,out);
            if(success==1){
                results[newPntCount]=out;
                newPntCount+=1;
                if(newPntCount==2) break;
            }
        }
        result1 = results[0];
        result2 = results[1];

        if(newPntCount==2) createIntersection(result1,result2,targetPrimNum);
    }
}
float bounds[]=primintrinsic(0,"bounds",@primnum);
int nearPrims[] = primfind(1,set(bounds[0],bounds[2],bounds[4]),set(bounds[1],bounds[3],bounds[5]));
foreach(int nearPrim;nearPrims){
    int check;
    vector pnt1;
    vector pnt2;
    getPrimIntersection(@primnum,nearPrim,check,pnt1,pnt2);
}

長っっっっっっ！:sweat:
これでコードによってオブジェクト同士の交線を求めることができました。
豚さん交線

まぁ、HIPファイルを添付するので見たい方はどうぞ。こちらです。

5
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
5
1