C# Vector3D.Cross方法代码示例

		/// <summary>
		/// Calculates the angle subtended by two line segments that meet at a vertex.
		/// </summary>
		/// <param name="start">The end of one of the line segments.</param>
		/// <param name="vertex">The vertex of the angle formed by the two line segments.</param>
		/// <param name="end">The end of the other line segment.</param>
		/// <returns>The angle subtended by the two line segments in degrees.</returns>
		public static double SubtendedAngle(PointF start, PointF vertex, PointF end)
		{
			Vector3D vertexPositionVector = new Vector3D(vertex.X, vertex.Y, 0);
			Vector3D a = new Vector3D(start.X, start.Y, 0) - vertexPositionVector;
			Vector3D b = new Vector3D(end.X, end.Y, 0) - vertexPositionVector;

			float dotProduct = a.Dot(b);

			Vector3D crossProduct = a.Cross(b);

			float magA = a.Magnitude;
			float magB = b.Magnitude;

			if (FloatComparer.AreEqual(magA, 0F) || FloatComparer.AreEqual(magB, 0F))
				return 0;

			double cosTheta = dotProduct/magA/magB;

			// Make sure cosTheta is within bounds so we don't
			// get any errors when we take the acos.
			if (cosTheta > 1.0f)
				cosTheta = 1.0f;

			if (cosTheta < -1.0f)
				cosTheta = -1.0f;

			double theta = Math.Acos(cosTheta)*(crossProduct.Z == 0 ? 1 : -Math.Sign(crossProduct.Z));
			double thetaInDegrees = theta/Math.PI*180;

			return thetaInDegrees;
		}

        public static double DistanceLineLine( Vector3D n1, Vector3D p1, Vector3D n2, Vector3D p2 )
        {
            // Check to make sure that neither of the normal vectors are degenerate.
            if( Tolerance.Zero( n1.MagSquared() ) || Tolerance.Zero( n2.MagSquared() ) )
                return double.NaN;

            Vector3D plane = n1.Cross( n2 );

            //	Case where the lines are parallel (magnitude of the cross product will be 0).
            if( Tolerance.Zero( plane.MagSquared() ) )
                return DistancePointLine( n1, p1, p2 );

            return DistancePointPlane( plane, p1, p2 );
        }

        /// <summary>
        /// Returns the two intersection points of two great circles.
        /// gc1 and gc2 are the great circle normal vectors.
        /// Fails if the input circles are the same.
        /// </summary>
        public static bool IntersectionGCGC( Vector3D gc1, Vector3D gc2, out Vector3D i1, out Vector3D i2 )
        {
            i1 = i2 = Vector3D.DneVector();

            // Direction vector of the intersection point of the two great circles.
            // NOTE there are actually two antipodal intersection points, +-I
            Vector3D I = gc1.Cross( gc2 );
            if( !I.Normalize() )
                return false;

            i1 = I;
            i2 = -I;
            return true;
        }

        /// <summary>
        /// NOTE: Not general, and assumes some things we know about this problem domain, 
        /// e.g. that c1 and c2 live on the same sphere of radius 1, and have two intersection points.
        /// </summary>
        public static void IntersectionCircleCircle( Vector3D sphereCenter, Circle3D c1, Circle3D c2, out Vector3D i1, out Vector3D i2 )
        {
            // Spherical analogue of our flat circle-circle intersection.
            // Spherical pythagorean theorem for sphere where r=1: cos(hypt) = cos(A)*cos(B)

            Circle3D clone1 = c1.Clone(), clone2 = c2.Clone();
            //clone1.Center -= sphereCenter;
            //clone2.Center -= sphereCenter;

            // Great circle (denoted by normal vector), and distance between the centers.
            Vector3D gc = clone2.Normal.Cross( clone1.Normal );
            double d = clone2.Normal.AngleTo( clone1.Normal );
            double r1 = clone1.Normal.AngleTo( clone1.PointOnCircle );
            double r2 = clone2.Normal.AngleTo( clone2.PointOnCircle );

            // Calculate distances we need.  So ugly!
            // http://www.wolframalpha.com/input/?i=cos%28r1%29%2Fcos%28r2%29+%3D+cos%28x%29%2Fcos%28d-x%29%2C+solve+for+x
            double t1 = Math.Pow( Math.Tan( d / 2 ), 2 );
            double t2 = Math.Cos( r1 ) / Math.Cos( r2 );
            double t3 = Math.Sqrt( (t1 + 1) * (t1 * t2 * t2 + 2 * t1 * t2 + t1 + t2 * t2 - 2 * t2 + 1) ) - 2 * t1 * t2;
            double x = 2 * Math.Atan( t3 / (t1 * t2 + t1 - t2 + 1) );
            double y = Math.Acos( Math.Cos( r1 ) / Math.Cos( x ) );

            i1 = clone1.Normal;
            i1.RotateAboutAxis( gc, x );
            i2 = i1;

            // Perpendicular to gc through i1.
            Vector3D gc2 = i1.Cross( gc );
            i1.RotateAboutAxis( gc2, y );
            i2.RotateAboutAxis( gc2, -y );
            i1 += sphereCenter;
            i2 += sphereCenter;

            /*
            // It would be nice to do the spherical analogue of circle-circle intersections, like here:
            // http://mathworld.wolfram.com/Circle-CircleIntersection.html
            // But I don't want to jump down that rabbit hole and am going to sacrifice some speed to use
            // my existing euclidean function.

            // Stereographic projection to the plane.  XXX - Crap, circles may become lines, and this isn't being handled well.
            Circle3D c1Plane = H3Models.BallToUHS( clone1 );
            Circle3D c2Plane = H3Models.BallToUHS( clone2 );
            if( 2 != Euclidean2D.IntersectionCircleCircle( c1Plane.ToFlatCircle(), c2Plane.ToFlatCircle(), out i1, out i2 ) )
                throw new System.Exception( "Expected two intersection points" );
            i1 = H3Models.UHSToBall( i1 ); i1 += sphereCenter;
            i2 = H3Models.UHSToBall( i2 ); i2 += sphereCenter;
            */
        }

        /// <summary>
        /// Given a curvature tensor, finds principal directions and curvatures.
        /// </summary>
        /// <param name="u">The u vector.</param>
        /// <param name="v">The v vector.</param>
        /// <param name="ku">ku.</param>
        /// <param name="kuv">kuv.</param>
        /// <param name="kv">kv.</param>
        /// <param name="normalNew">The new normal.</param>
        /// <param name="pDirMax">The maximum principle curvature direction.</param>
        /// <param name="pDirMin">The minimum principle curvature direction.</param>
        /// <param name="kMax">The maximum principle curvature.</param>
        /// <param name="kMin">The minimum principle curvature.</param>
        public static void DiagonalizeCurvature(Vector3D u, Vector3D v, double ku, double kuv, double kv, Vector3D normalNew,
            out Vector3D pDirMax, out Vector3D pDirMin, out double kMax, out double kMin)
        {
            Vector3D uRotated, vRotated;

            RotateCoordinateSystem(u, v, normalNew, out uRotated, out vRotated);

            double c = 1.0 ;
            double s = 0.0 ;
            double t = 0.0 ;

            // Jacobi rotation to diagonalize
            if (kuv != 0.0f)
            {
                double h = 0.5f * (kv - ku) / kuv;

                if (h < 0.0f)
                {
                    t = 1.0f / (h - (float)Math.Sqrt(1.0f + h * h));
                }
                else
                {
                    t = 1.0f / (h + (float)Math.Sqrt(1.0f + h * h));
                }

                c = 1.0f / (float)Math.Sqrt(1.0f + t * t);
                s = t * c;
            }

            kMax = ku - t * kuv;
            kMin = kv + t * kuv;

            if (Math.Abs(kMax) >= Math.Abs(kMin))
            {
                pDirMax = c * uRotated - s * vRotated;
            }
            else
            {
                double temp = kMin;
                kMin = kMax;
                kMax = temp;
                pDirMax = s * uRotated + c * vRotated;
            }

            pDirMin =  normalNew.Cross(pDirMax);
        }

        // 内心
        public void ComputeIncircleCenter(out Vector3D center, out double radius)
        {
            center = Vector3D.Zero;
            radius = 0;
            double a, b, c;

            c = ComputeEdgeLength(1);
            b = ComputeEdgeLength(3);
            a = ComputeEdgeLength(2);

            center.x = (a * this.A.x + b * this.B.x + c * this.C.x) / (a + b + c);
            center.y = (a * this.A.y + b * this.B.y + c * this.C.y) / (a + b + c);
            center.z = (a * this.A.z + b * this.B.z + c * this.C.z) / (a + b + c);

            radius = center.Cross(C - A).Length();

        }

		public void TestCross()
		{
			Vector3D v1 = new Vector3D(2.2F, -6.1F, 7.4F);
			Vector3D v2 = new Vector3D(-3.8F, 3.7F, 4.1F);
			Vector3D result = new Vector3D(-52.39F, -37.14F, -15.04F);

			Assert.IsTrue(Vector3D.AreEqual(v1.Cross(v2), result));
		}

 public static Circle3D OrthogonalCircle( Vector3D v1, Vector3D v2 )
 {
     Vector3D center;
     double rad;
     OrthogonalCircle( v1, v2, out center, out rad );
     Vector3D normal = v1.Cross( v2 );
     return new Circle3D { Center = center, Normal = normal, Radius = rad };
 }

			private static Matrix ImageOrientationPatientToMatrix(ImageOrientationPatient orientation)
			{
				Vector3D xOrient = new Vector3D((float) orientation.RowX, (float) orientation.RowY, (float) orientation.RowZ);
				Vector3D yOrient = new Vector3D((float) orientation.ColumnX, (float) orientation.ColumnY, (float) orientation.ColumnZ);
				Vector3D zOrient = xOrient.Cross(yOrient);

				Matrix orientationMatrix = Math3D.OrientationMatrixFromVectors(xOrient, yOrient, zOrient);
				return orientationMatrix;
			}

        /// <summary>
        /// Projects a curvature tensor from an old basis to a new one.
        /// </summary>
        /// <param name="uOld">Old u.</param>
        /// <param name="vOld">Old v.</param>
        /// <param name="kuOld">Old ku.</param>
        /// <param name="kuvOld">Old kuv.</param>
        /// <param name="kvOld">Old kv.</param>
        /// <param name="uNew">New u</param>
        /// <param name="vNew">New v.</param>
        /// <param name="kuNew">New ku.</param>
        /// <param name="kuvNew">New kuv.</param>
        /// <param name="kvNew">New kv.</param>
        public static void ProjectCurvature(Vector3D uOld, Vector3D vOld, float kuOld, float kuvOld, float kvOld, Vector3D uNew, Vector3D vNew,
            out float kuNew, out float kuvNew, out float kvNew)
        {
            Vector3D uNewRotated, vNewRotated;

            RotateCoordinateSystem(uNew, vNew,  uOld.Cross(vOld), out uNewRotated, out vNewRotated);

            float u1 =(float) uNewRotated.Dot(uOld);
            float v1 = (float)uNewRotated.Dot(vOld);
            float u2 = (float)vNewRotated.Dot(uOld);
            float v2 = (float)vNewRotated.Dot(vOld);

            kuNew  = kuOld * u1 * u1 + kuvOld * (2.0f    * u1 * v1) + kvOld * v1 * v1;
            kuvNew = kuOld * u1 * u2 + kuvOld * (u1 * v2 + u2 * v1) + kvOld * v1 * v2;
            kvNew  = kuOld * u2 * u2 + kuvOld * (2.0f    * u2 * v2) + kvOld * v2 * v2;
        }

 /// <summary>
 /// Calculate the sine of the angle between two vectors.
 /// </summary>
 private static double SinAngle( Vector3D p1, Vector3D p2 )
 {
     double sinA = (p1.Cross( p2 )).Abs() / (p1.Abs() * p2.Abs());
     return Clamp( sinA );
 }

 /// <summary>
 /// Return the normal of the great circle defined by the two vectors. 
 /// Return false if we fail to get the normal (happens if p1 = p2).
 /// </summary>
 static bool GetGCNormal( Vector3D p1, Vector3D p2, out Vector3D normal )
 {
     normal = p1.Cross( p2 );
     if( !normal.Normalize() )
     {
         normal = Vector3D.DneVector();
         return false;
     }
     return true;
 }

        public static IPresentationImage PresentationImageFromPositionOrientation(
            ImagePositionPatient imagePositionPatient,
            ImageOrientationPatient imageOrientationPatient,
            IDisplaySet displaySet,
            string frameOfReferenceUid)
        {
            var point = new Vector3D(
                (float) imagePositionPatient.X,
                (float) imagePositionPatient.Y,
                (float) imagePositionPatient.Z);

            if (displaySet != null)
            {
                var firstSop = displaySet.PresentationImages.OfType<IImageSopProvider>().FirstOrDefault();
                // Match Frame of Reference UID, if present
                if (firstSop == null || string.IsNullOrEmpty(frameOfReferenceUid) ||
                    string.IsNullOrEmpty(firstSop.Frame.FrameOfReferenceUid) || frameOfReferenceUid == firstSop.Frame.FrameOfReferenceUid)
                {
                    foreach (IPresentationImage image in displaySet.PresentationImages)
                    {
                        var sop = image as IImageSopProvider;
                        if (sop != null)
                        {
                            Vector3D planeRow = new Vector3D(
                                (float) sop.Frame.ImageOrientationPatient.RowX,
                                (float) sop.Frame.ImageOrientationPatient.RowY,
                                (float) sop.Frame.ImageOrientationPatient.RowZ).Normalize();
                            Vector3D planeColumn = new Vector3D(
                                (float) sop.Frame.ImageOrientationPatient.ColumnX,
                                (float) sop.Frame.ImageOrientationPatient.ColumnY,
                                (float) sop.Frame.ImageOrientationPatient.ColumnZ).Normalize();
                            Vector3D planeNormal = planeRow.Cross(planeColumn);

                            var planePosition = new Vector3D(
                                (float) sop.Frame.ImagePositionPatient.X,
                                (float) sop.Frame.ImagePositionPatient.Y,
                                (float) sop.Frame.ImagePositionPatient.Z);
                            if (PointIsInPlane(point, planeNormal, planePosition))
                                return image;
                        }
                    }
                }
            }

            return null;
        }

        public IList<Vector3D> ComputeCurve(
            Vector3D start, 
            Vector3D stop, 
            double granularity)
        {
            if (granularity <= 0.0)
            {
                throw new ArgumentOutOfRangeException("granularity", "Granularity must be greater than zero.");
            }

            Vector3D normal = start.Cross(stop).Normalize();
            double theta = start.AngleBetween(stop);
            int n = Math.Max((int)(theta / granularity) - 1, 0);
            
            List<Vector3D> positions = new List<Vector3D>(2 + n);

            positions.Add(start);

            for (int i = 1; i <= n; ++i)
            {
                double phi = (i * granularity);

                positions.Add(ScaleToGeocentricSurface(start.RotateAroundAxis(normal, phi)));
            }

            positions.Add(stop);

            return positions;
        }

        /// <summary>
        /// Used to get the perpendicular to the polyline at a point p2, given adjacent points.
        /// </summary>
        private static Vector3D CurvePerp( Vector3D p1, Vector3D p2, Vector3D p3 )
        {
            // Just use p1 and p3 for now, close to the 3 point rule.
            // http://math.depaul.edu/mash/optnum.pdf

            Vector3D perpendicular = p1.Cross( p3 );
            //Vector3D perpendicular = axis.Cross( new Vector3D( 0, 0, 1 ) );
            if( !perpendicular.Normalize() )
            {
                // This can happen if p1 and p3 are collinear with origin.
                Vector3D axis = p1 - p3;
                perpendicular = axis.Cross( new Vector3D( 0, 1, 0 ) );
                perpendicular.Normalize();
            }

            return perpendicular;
        }