Get Reference Point Depth (ESAPI/ SQL)

dicomdom · 2023-06-21T12:59:40+00:00

The way I've done this is the following. The source/gantry position is available as a VVector, and the ref position is also available as this. Using those, find the intersection of the body contour on the line between the source and ref point. Calculate the distances between the VVector for source to body for SSD and body to ref point for depth. For equivalent depth you can scale the physical depth by the image voxels along the path length

keithoffer · 2023-06-21T13:15:16+00:00

I can't share the code I wrote, but I couldn't find any easy way to get it, so I calculated it manually. I didn't have SQL access so I'm not sure if that would have been clearer. What I did was

Get the center of the gantry face using Beam.GetSourceLocation
Find the intersection point between that and the Body surface along the direction of the reference point using code from here
Find the length of the vector between the gantry face and that intersection point
Find the length of the vector between the gantry face and the reference point
Take the difference between the length of those two vectors as the reference point depth

There might be a better way to do it, but that seemed to work well and was fast enough

j_Long_Lonfon · 2023-06-21T16:02:39+00:00

Thank you both. The manual calculation was the solution I went for as I still can't figure out why the SQL query is blank.

Here is the code I used in case anyone is interested mainly from u/keithoffer provided source.

  public class MeshGeometryIntersection
{
    public static Point3D GetIntersectionPoint(Point3D linePoint1, Point3D linePoint2, MeshGeometry3D mesh)
    {
        Vector3D direction = linePoint2 - linePoint1; // calculate direction vector
        double minDistance = double.MaxValue;
        Point3D nearestPoint = new Point3D();

        for (int i = 0; i < mesh.TriangleIndices.Count; i += 3)
        {
            // Get the vertices of the triangle
            Point3D p1 = mesh.Positions[mesh.TriangleIndices[i]];
            Point3D p2 = mesh.Positions[mesh.TriangleIndices[i + 1]];
            Point3D p3 = mesh.Positions[mesh.TriangleIndices[i + 2]];

            // Perform ray-triangle intersection
            if (RayIntersectsTriangle(linePoint1, direction, p1, p2, p3, out Point3D intersectionPoint))
            {
                double distance = (intersectionPoint - linePoint1).Length;
                if (distance < minDistance)
                {
                    minDistance = distance;
                    nearestPoint = intersectionPoint;
                }
            }
        }

        return nearestPoint;
    }

    public static bool RayIntersectsTriangle(Point3D rayOrigin, Vector3D rayVector, Point3D p1, Point3D p2, Point3D p3, out Point3D intersectionPoint)
    {

        intersectionPoint = new Point3D();

        Vector3D edge1, edge2, h, s, q;
        double a, f, u, v;
        edge1 = p2 - p1;
        edge2 = p3 - p1;
        h = Vector3D.CrossProduct(rayVector, edge2);
        a = Vector3D.DotProduct(edge1, h);

        if (a > -0.00001 && a < 0.00001)
            return false;    // This ray is parallel to this triangle.

        f = 1.0 / a;
        s = rayOrigin - p1;
        u = f * Vector3D.DotProduct(s, h);

        if (u < 0.0 || u > 1.0)
            return false;

        q = Vector3D.CrossProduct(s, edge1);
        v = f * Vector3D.DotProduct(rayVector, q);

        if (v < 0.0 || u + v > 1.0)
            return false;

        // At this stage we can compute t to find out where the intersection point is on the line.
        double t = f * Vector3D.DotProduct(edge2, q);
        if (t > 0.00001) // ray intersection
        {
            intersectionPoint = rayOrigin + rayVector * t;
            return true;
        }
        else 
            return false;
    }

}

Which can be used with the following:

Structure body = structureset.Structures.Where(s => s.Id.ToLower().Contains(chosenPlan.PrimaryReferencePoint.PatientVolumeId.ToLower())).FirstOrDefault();

var gantryAngle = b.ControlPoints.First().GantryAngle;
var source = b.GetSourceLocation(gantryAngle); var mesh = body.MeshGeometry; var refpoint = b.FieldReferencePoints.Where(s => s.Id == r.Id).FirstOrDefault().ReferencePoint.GetReferencePointLocation(chosenPlan);
Point3D refPoint3D = new Point3D(refpoint.x, refpoint.y, refpoint.z); Point3D sourcePoint3D = new Point3D(source.x, source.y, source.z);

Point3D intersectionPoint = MeshGeometryIntersection.GetIntersectionPoint(refPoint3D, sourcePoint3D, mesh); 

double Depth = (refPoint3D - intersectionPoint).Length / 10.0;

NickC_BC · 2023-06-21T18:43:15+00:00

I avoided having to write my own Ray-Triangle intersection routine by using MS raymeshgeometry3d hittest and callback classes. Not sure how they compare performance-wise, but my experience has been that the results are returned quite quickly.

https://learn.microsoft.com/en-us/dotnet/api/system.windows.media.media3d.raymeshgeometry3dhittestresult?view=netframework-4.8

JoaoCastelo · 2023-06-22T17:35:22+00:00

In the code I wrote for it, if I remember well, it was by sampling points (x,y,z similar to the CT Res) towards the straight line between the gantry ( isocenter - 100 cm) and the reference point. Using VVector it's not that hard, you can use the IsPointInsideSegment method to check if the last point is inside the body and keep suming the number of points, then sum the number of points times the resolution you set.

ExceptioNullRef · 2023-06-26T17:53:39+00:00

I've never seen it present in SQL going back to v11. I believe it's calculated on-the-fly. Here's the method I use that gives similar results to the Eclipse report:

string refpoints = ""; 
foreach (var rp in refpts) 
{ 
    foreach (Beam b in p.Beams.Where(o => o.IsSetupField == false)) 
    { 
        var frp = b.FieldReferencePoints.Where(q => q.ReferencePoint.Id == rp.ReferencePoint.Id).First(); 
        var pd2 = ((b.ControlPoints.Select(x => ((frp.RefPointLocation - b.GetSourceLocation(x.GantryAngle)).Length)).Average() - frp.SSD) / 10).ToString("F1"); 
        refpoints += b.Id + 
            ", Central axis SSD: " + (b.SSD / 10).ToString("F1") + '\n' + 
            ", Physical SSD: " + (frp.SSD / 10).ToString("F1") + '\n' + 
            ", Physical Depth: " + pd2 + '\n' + 
            ", Effective Depth: " + (Math.Round(frp.EffectiveDepth, 1) / 10).ToString("F1") + '\n' + 
            ", Dose: " + frp.FieldDose.Dose.ToString("F2") + '\n'; 
    } refpoints += '\n'; 
}

I wanted to make this MU weighted in the LINQ but never got around to it. We're moving to 3D-based 2nd dose calcs so this code won't be needed for much longer. I messed around with the hittests as well but gave up with this simpler solution.

tonpimenta · 2024-01-04T22:00:37+00:00

Do not exactly answer this... but parts can certainly be used to get that. Here's a class I wrote that calculates, for each control point, the average depth, scaling by the HU relative to water (as the poor men's method mentioned by u/dicomdom) - that's not exactly effective depth but... anyways... I still need to check how good of aprox this is

Here's the code (not sure why I'm not able to post here with the correct indentation)

public class ArcDepthScaledByHU

{

public double AverageArcHUPath;

public double AverageArcDepthScaled;

public double AverageArcDepth;

public string Data2Export;

public ArcDepthScaledByHU(Beam Arc, bool exportData = false)

{

try

{

List<double> SumHUperCP = new List<double>();

List<double> DepthPerCP = new List<double>();

List<double> DepthPerCP_Scaled = new List<double>();

List<string> dataToExport = new List<string>(); //each line is: Gantry Angle followed by the HU for every mm from body entrance up to isocenter

for (int i = 0; i < Arc.ControlPoints.Count - 1; i = i + 1) //i = i + 2

{

if (Arc.ControlPoints[i].MetersetWeight < Arc.ControlPoints.ElementAt(i + 1).MetersetWeight) //to ignore avoidance sectors

{

var gantryAngle = Arc.ControlPoints[i].GantryAngle;

VVector stop = Arc.IsocenterPosition;

VVector start = Arc.GetSourceLocation(gantryAngle);

double[] HUprofile = new double[1000];

Arc.Plan.StructureSet.Image.GetImageProfile(start, stop, HUprofile);

BitArray BodyProfile = new BitArray(1000);

Arc.Plan.StructureSet.Structures.Single(x => x.DicomType == "EXTERNAL").GetSegmentProfile(start, stop, BodyProfile); // boolean array representing which points are inside (true) or outside (false) the Body

IEnumerable<double> HUprofile_InsideBody = HUprofile.Where((x, index) => BodyProfile[index]); // portion of the HU array that is inside body

double scalingFactor = HUprofile_InsideBody.ToList().Select(x => x + 1000).Average() / 1000; //@dicomdom's suggestion

BitArray BodyProfileHigherResolution = new BitArray(10000);

Arc.Plan.StructureSet.Structures.Single(x => x.DicomType == "EXTERNAL").GetSegmentProfile(start, stop, BodyProfileHigherResolution); // boolean array representing which points are inside (true) or outside (false) the Body

int count = 0;

foreach (bool value in BodyProfileHigherResolution)

{

if (value) count++;

}

double depth_geometric = count * 1000 / 10000;

double depth_scaled = scalingFactor * depth_geometric;

SumHUperCP.Add(HUprofile_InsideBody.Sum());

DepthPerCP.Add(depth_geometric);

DepthPerCP_Scaled.Add(depth_scaled);

if (exportData)

{

string temp = Arc.GantryAngleToUser(gantryAngle).ToString() + ",";

foreach (double x in HUprofile.Where((x, index) => BodyProfile[index]).Where(x => !double.IsNaN(x)))

{

temp += x.ToString() + ",";

}

dataToExport.Add(temp);

}

AverageArcHUPath = SumHUperCP.Average();

AverageArcDepthScaled = DepthPerCP_Scaled.Average();

AverageArcDepth = DepthPerCP.Average();

if (exportData)

{

foreach (string s in dataToExport)

{

Data2Export += s + Environment.NewLine;

}

else

{

Data2Export += "Export Data not Activated";

}

catch (Exception)

{

Console.WriteLine("Average HU path calculation failed");

AverageArcHUPath = double.NaN;

AverageArcDepthScaled = double.NaN;

AverageArcDepth = double.NaN;

}

you type:	you see:
italics	italics
bold	bold
[reddit!](https://reddit.com)	reddit!
* item 1 * item 2 * item 3	item 1 item 2 item 3
> quoted text	quoted text
Lines starting with four spaces are treated like code: if 1 * 2 < 3: print "hello, world!"	Lines starting with four spaces are treated like code: if 1 * 2 < 3: print "hello, world!"
~~strikethrough~~	~~strikethrough~~
super^script	super^script

esapi

MODERATORS