<<>>
<<ImplicitTubeImpl.mesa>>
<<Copyright Š 1985 by Xerox Corporation.  All rights reserved.>>
<<Bloomenthal, September 6, 1992 12:01 pm PDT>>

DIRECTORY G3dBasic, G3dDraw, G3dFunction, G3dMatrix, G3dSpline, G3dTube, G3dVector, Imager, ImplicitDefs, ImplicitTube, RealFns, 
ImplicitSurface;

ImplicitTubeImpl: CEDAR MONITOR
IMPORTS G3dDraw, G3dFunction, G3dSpline, G3dTube, G3dVector, RealFns, ImplicitSurface
EXPORTS ImplicitTube

~ BEGIN

Pair:                TYPE ~ G3dBasic.Pair;
Triple:            TYPE ~ G3dBasic.Triple;
Matrix:            TYPE ~ G3dMatrix.Matrix;
Tube:                TYPE ~ G3dTube.Tube;
TubeProc:        TYPE ~ G3dTube.TubeProc;
XSection:            TYPE ~ G3dTube.XSection;
NearSegment:    TYPE ~ G3dVector.NearSegment;
Context:            TYPE ~ Imager.Context;
Sample:            TYPE ~ ImplicitDefs.Sample;
DistanceMode:    TYPE ~ ImplicitDefs.DistanceMode;
ValueProc:        TYPE ~ ImplicitDefs.ValueProc;

PI:                    REAL ~ 3.1415926535;

PrepareTube: PUBLIC PROC [tube: Tube, distanceMode: DistanceMode, measureMode: ATOM]
~ {
IF NOT tube.xSectionsValid THEN G3dTube.MakeXSections[tube, FALSE, TRUE];
IF NOT tube.lengthsValid THEN G3dTube.SetLengths[tube];
IF NOT tube.planeValid THEN G3dTube.SetPlanes[tube];
SELECT measureMode FROM
$Sum =>
G3dTube.SetDescendantRadii[tube, IF distanceMode=inverse THEN linear ELSE square];
ENDCASE;
};

StartPoint: PUBLIC PROC [
tube: Tube,
valueProc: ValueProc,
threshold: REAL Ź 1.0,
measureMode: ATOM Ź $Unspecified,
spread: REAL Ź 0.75,
clientData: REF ANY Ź NIL]
RETURNS [tt: Triple]
~ {
t: REAL ~ 0.5;
pIn, pOut, v: Triple;
factor: REAL;
factor Ź SELECT measureMode FROM
$Segment => 2.0/spread,
--($Sum)-- ENDCASE => 1.5;    -- 1.5 instead of 1.0, for some fudge
pIn Ź G3dSpline.Position[tube.spline, t];
v Ź G3dSpline.Velocity[tube.spline, t];
v Ź G3dVector.Unit[G3dVector.Ortho[v]];
pOut Ź G3dVector.ScaleRay[[pIn, v], factor*G3dTube.Radius[tube, t]];
tt Ź ImplicitSurface.SurfacePoint[pIn, pOut, valueProc, threshold, clientData];
};

SampleTube: PUBLIC PROC [
q: Triple,
tube: Tube,
measureMode: ATOM,
distanceMode: DistanceMode,
tolerance, spread: REAL]
RETURNS [sample: Sample]
~ {
s: Sample Ź sample Ź
SampleTubeBranch[q, tube, measureMode, distanceMode, tolerance, spread];
nearest: Tube Ź G3dTube.NearestTube[q, tube].tube;
IF nearest.nearSpline.distance # 0.0
THEN {                                                -- useful sample information:
v: Triple ~ G3dVector.Sub[q, nearest.nearSpline.point];
radius: REAL ~ G3dTube.Radius[nearest, nearest.nearSpline.t];
d: REAL ~ nearest.nearSpline.distance;
sample.near Ź G3dVector.ScaleRay[[nearest.nearSpline.point, v], radius/d];
}
ELSE sample.near Ź nearest.nearSpline.point;
};

SampleTubeBranch: PUBLIC PROC [
q: Triple,
tube: Tube,
measureMode: ATOM,
distanceMode: DistanceMode,
tolerance, spread: REAL]
RETURNS [Sample]
~ {
segment: Sample Ź
SampleTubeSegment[q, tube, measureMode, distanceMode, tolerance, spread];
branches: Sample Ź
SampleTubeBranches[q, tube, measureMode, distanceMode, tolerance, spread];
<<Note about the following line: if tube has no next and no branches and q is too far away from tube, then both segment.value and 
branches.value may be zero; in this case, segment.near, segment.vector, etc. are defined, but branches.near, etc. are not; thus the 
benefit of the inequality goes to segment.>>
RETURN[IF branches.value > segment.value THEN branches ELSE segment];
};

<<SampleTubeBranch: PUBLIC PROC [>>
<<q: Triple,>>
<<tube: Tube,>>
<<sampleMode: SampleMode]>>
<<RETURNS [Sample]>>
<<~ {>>
<<If we want to debug the contribution of branches only:>>
<<RETURN[IF tube.prev = NIL THEN SampleTubeBranches[q, tube, sampleMode] ELSE SampleTubeSegment[q, tube, sampleMode]];>>
<<};>>
<<>>
SampleTubeBranches: PUBLIC PROC [
q: Triple,
tube: Tube,
measureMode: ATOM,
distanceMode: DistanceMode,
tolerance, spread: REAL]
RETURNS [sample: Sample]
~ {
tubeProc: TubeProc ~ {
s: Sample Ź
SampleTubeBranch[q, tube, measureMode, distanceMode, tolerance, spread];
IF s.value > sample.value THEN sample Ź s;
samples[nBranches] Ź s;
nBranches Ź nBranches+1;
};
nBranches: NAT Ź 0;
samples: ARRAY [0..50) OF Sample;
G3dTube.ApplyToBranches[tube, tubeProc];
sample.value Ź 0.0;
SELECT measureMode FROM
$Segment => {
FOR n: NAT IN [0..nBranches) DO
s: Sample Ź samples[n];
FOR nn: NAT IN (n..nBranches) DO
ss: Sample Ź samples[nn];
IF s.near = ss.near
THEN sample.value Ź sample.value+s.value
ELSE {
aveRadius: REAL Ź 0.5*(s.radius+ss.radius);
aveDistance: REAL Ź 0.5*(s.distance+ss.distance);
aveRatio: REAL Ź Ratio[aveDistance, aveRadius, spread];
SELECT aveRatio FROM
< 2.0 => {
near: NearSegment Ź
G3dVector.NearestToSegment[s.near, ss.near, q];
lineDistance: REAL Ź SELECT near.w0 FROM
0.0 => ss.distance,
1.0 => s.distance,
ENDCASE    => G3dVector.Distance[near.point, q];
lineValue: REAL Ź
G3dFunction.Wyvill[Ratio[lineDistance, aveRadius, spread]];
IF aveRatio < 1.0
THEN sample.value Ź sample.value+lineValue
ELSE {
f: REAL Ź 2.0-aveRatio;
sample.value Ź sample.value+f*lineValue+(1.-f)*(s.value+ss.value);
};
};
ENDCASE => sample.value Ź sample.value+s.value+ss.value;
};
ENDLOOP;
ENDLOOP;
IF nBranches > 2 THEN sample.value Ź 2.0*sample.value/(nBranches*(nBranches-1));
};
--($Sum)-- ENDCASE =>
FOR n: NAT IN [0..nBranches) DO
sample.value Ź sample.value+samples[n].value;
ENDLOOP;
};

SampleTubeSegment: PUBLIC PROC [
q: Triple,
tube: Tube,
measureMode: ATOM,
distanceMode: DistanceMode,
tolerance, spread: REAL]
RETURNS [sample: Sample]
~ {
SetBasics: PROC ~ {
tube.nearSpline Ź G3dSpline.NearestPoint[q, tube.spline, 0.0, 1.0, tolerance];
sample.near Ź tube.nearSpline.point;
sample.vector Ź G3dVector.Sub[sample.near, q];
sample.radius Ź G3dTube.Radius[tube, tube.nearSpline.t];
sample.squareDistance Ź G3dVector.Square[sample.vector];
};
SetDistance: PROC ~ {
sample.distance Ź RealFns.SqRt[sample.squareDistance];
sample.distanceSet Ź TRUE;
};
sample.refAny Ź tube;
SELECT measureMode FROM
$Segment => {
SetBasics[];
SetDistance[];
sample.ratio Ź Ratio[sample.distance, sample.radius, spread];
sample.value Ź G3dFunction.Wyvill[sample.ratio];
};
--($Sum)-- ENDCASE => {
SetBasics[];
IF distanceMode = inverse THEN SetDistance[];
sample.value Ź IF distanceMode = inverse
THEN sample.radius/MAX[0.0001, sample.distance]
ELSE sample.radius*sample.radius/MAX[0.0001, sample.squareDistance];
};
};

Ratio: PROC [distance, radius, spread: REAL] RETURNS [REAL] ~ {
RETURN[IF radius = 0.0 THEN 2.0 ELSE spread*(distance-radius)/radius];
};

TextureOfTube: PUBLIC PROC [q: Triple, tube: Tube] RETURNS [texture: Pair] ~ {
nearest: Tube Ź G3dTube.NearestTube[q, tube, FALSE].tube;
xSection: XSection Ź G3dTube.GetXSection[nearest, nearest.nearSpline.t];
vector: Triple Ź G3dVector.Sub[q, xSection.frame.position];
xProj: Triple Ź G3dVector.Project[vector, xSection.frame.triad.n];
yProj: Triple Ź G3dVector.Project[vector, xSection.frame.triad.b];
proj: Triple Ź G3dVector.Unit[G3dVector.Add[xProj, yProj]];
texture.x Ź G3dVector.AngleBetween[proj, xSection.frame.triad.n, FALSE]/PI;
texture.y Ź SELECT nearest.nearSpline.t FROM
0.0 => tube.length0-
G3dVector.Length[G3dVector.Project[vector, xSection.frame.triad.v]],
1.0 => tube.length1+
G3dVector.Length[G3dVector.Project[vector, xSection.frame.triad.v]],
ENDCASE => nearest.length0+xSection.length;
};

DiagramTube: PUBLIC PROC [
context: Context,
view: Matrix,
viewport: G3dMatrix.Viewport,
q: Triple,
tube: Tube,
measureMode: ATOM,
distanceMode: DistanceMode,
tolerance, spread: REAL]
~ {
Info: TYPE ~ RECORD [tube: Tube, branches: BOOL, sample: Sample];
InfoTubeBranch: PROC [q: Triple, tube: Tube] RETURNS [Info] ~ {
segment: Info Ź InfoTubeSegment[q, tube];
branches: Info Ź InfoTubeBranches[q, tube];
RETURN[IF branches.sample.value > segment.sample.value THEN branches ELSE segment];
};
InfoTubeBranches: PROC [q: Triple, tube: Tube] RETURNS [info: Info] ~ {
s: Sample ŹSampleTubeSegment[q, tube, measureMode, distanceMode, tolerance, spread];
RETURN[[tube, TRUE, s]];
};
InfoTubeSegment: PROC [q: Triple, tube: Tube] RETURNS [Info] ~ {
s: Sample ŹSampleTubeSegment[q, tube, measureMode, distanceMode, tolerance, spread];
RETURN[[tube, FALSE, s]];
};
IF tube # NIL THEN SELECT measureMode FROM
$Segment => {
info: Info Ź InfoTubeBranch[q, tube];
IF info.branches
THEN {
nBranches: NAT Ź G3dTube.NBranches[info.tube]+1;
FOR n: NAT IN [0..nBranches) DO
b0: Tube Ź G3dTube.GetBranch[info.tube, n];
p0: Triple Ź
G3dSpline.NearestPoint[q, b0.spline, 0.0, 1.0, tolerance].point;
FOR nn: NAT IN (n..nBranches) DO
b1: Tube Ź G3dTube.GetBranch[info.tube, nn];
p1: Triple Ź
G3dSpline.NearestPoint[q, b1.spline, 0.0, 1.0, tolerance].point;
x: Triple Ź G3dVector.NearestToSegment[p0, p1, q].point;
G3dDraw.Segment[context, p0, p1, view, viewport, solid];
G3dDraw.Segment[context, q, x, view, viewport, solid];
G3dDraw.Segment[context, q, p0, view, viewport, dotted];
G3dDraw.Segment[context, q, p1, view, viewport, dotted];
ENDLOOP;
ENDLOOP;
}
ELSE {
x: Triple Ź G3dTube.NearestTube[q, info.tube].tube.nearSpline.point;
G3dDraw.Segment[context, q, x, view, viewport, solid];
};
};
ENDCASE;
};

END.
..

TextureOfTube: PUBLIC PROC [q: Triple, tube: Tube] RETURNS [texture: Pair] ~ {
Info: TYPE ~ RECORD [inverseDistance: REAL, texture: Pair];
tubeProc: TubeProc ~ {
xSection: XSection Ź G3dTube.GetXSection[tube, tube.nearSpline.t];
vector: Triple Ź G3dVector.Sub[q, frame.position];
xProj: Triple Ź G3dVector.Project[vector, xSection.frame.triad.n];
yProj: Triple Ź G3dVector.Project[vector, xSection.frame.triad.b];
proj: Triple Ź G3dVector.Unit[G3dVector.Add[xProj, yProj]];
tx: REAL Ź G3dVector.AngleBetween[proj, xSection.frame.triad.n, FALSE]/PI;
ty: REAL Ź SELECT tube.nearSpline.t FROM
0.0 => tube.length0-
G3dVector.Length[G3dVector.Project[vector, xSection.frame.triad.v]],
1.0 => tube.length1+
G3dVector.Length[G3dVector.Project[vector, xSection.frame.triad.v]],
ENDCASE => tube.length0+frame.length;
info[nBranches] Ź [1.0/tube.nearSpline.distance, [tx, ty]];
nBranches Ź nBranches+1;
};
nBranches: NAT Ź 0;
info: ARRAY [0..50) OF Info;
nearest: Tube Ź G3dTube.NearestTube[q, tube, FALSE];
G3dTube.ApplyToSiblings[nearest, tubeProc];
IF nBranches = 1
THEN RETURN[info[0].texture]
ELSE {
inverseDistanceSum: REAL Ź 0.0;
texture Ź [0.0, 0.0];
FOR n: NAT IN [0..nBranches) DO
inverseDistanceSum Ź inverseDistanceSum+info[n].inverseDistance;
ENDLOOP;
FOR n: NAT IN [0..nBranches) DO
texture Ź Vector2.Add[texture, Vector2.Mul[info[n].texture, info[n].inverseDistance]];
ENDLOOP;
texture Ź Vector2.Div[texture, inverseDistanceSum];
};
};

ValueOfTubeBranches: PUBLIC PROC [
q: Triple,
tube: Tube,
distanceMode: DistanceMode,
tolerance: REAL Ź 0.0]
RETURNS [sample: Sample]
~ {
nBranches: NAT Ź 0;
vSum: Triple Ź [0.0, 0.0, 0.0];
maxValue: REAL Ź -100000.0;
vectors: ARRAY [0..50) OF Triple;
samples: ARRAY [0..50) OF Sample;
weight: BOOL Ź distanceMode = wtInverse OR distanceMode = wtInvSq;
tubeProc: G3dTube.TubeProc ~ {
IF tube # NIL THEN {
s: Sample Ź ValueOfTubeBranch[q, tube, distanceMode, tolerance];
IF s.value > maxValue THEN {maxValue Ź s.value; sample Ź s};
IF (distanceMode=wtInverse OR distanceMode=wtInvSq) AND NOT s.distanceSet THEN {
s.distance Ź Real.SqRt[s.squareDistance];
s.distanceSet Ź TRUE;
};
samples[nBranches] Ź s;
nBranches Ź nBranches+1;
};
};
G3dTube.ApplyToBranches[tube, tubeProc];
SELECT distanceMode FROM
wtInverse, wtInvSq => {
<<Now a miracle occurs.>>
FOR n: NAT IN [0..nBranches) DO
v: Triple Ź G3dVector.Sub[samples[n].nearest, q];
vectors[n] Ź G3dVector.Div[v, samples[n].distance];
vSum Ź G3dVector.Add[vSum, vectors[n]];
ENDLOOP;
vSum Ź G3dVector.Unit[vSum];
FOR n: NAT IN [0..nBranches) DO
f: REAL ~ 0.5*(1.0-G3dVector.Dot[vSum, vectors[n]]);
sample.value Ź sample.value+f*samples[n].value;
ENDLOOP;
};
ENDCASE => {
FOR n: NAT IN [0..nBranches) DO
sample.value Ź sample.value+samples[n].value;
ENDLOOP;
};
};

ValueOfTubeSegment: PUBLIC PROC [
q: Triple,
tube: Tube,
distanceMode: DistanceMode,
tolerance: REAL Ź 0.0]
RETURNS [s: Sample]
~ {
radius: REAL ~ G3dTube.Radius[tube, tube.nearSpline.t];
s Ź ImplicitValue.ValueOfPoint[q, tube.nearSpline.point, radius, distanceMode, threshold];
};

mode: SampleMode.integral => {
IF tube.length = 0.0 THEN G3dTube.SetSectionLengths[tube];
{
Val: PROC RETURNS [REAL] ~ {
distance: REAL Ź G3dVector.Length[[diff[0][0], diff[0][1], diff[0][2]]];
RETURN[Strength[Ratio[distance, radius, mode.spread]]];
};
radius: REAL Ź tube.r0;
diff: Matrix Ź G3dMatrix.ObtainMatrix[];
nSteps: NAT Ź Real.Round[tube.length/mode.stepSize];
dRadius: REAL Ź (tube.r1-tube.r0)/REAL[MAX[1, nSteps]];
diff^ Ź tube.spline^;
diff[3] Ź [diff[3][0]-q.x, diff[3][1]-q.y, diff[3][2]-q.z, 0.0];
diff Ź G3dSpline.ForwardDifference[diff, nSteps];
sample.value Ź Val[];
FOR n: NAT IN [0..nSteps) DO
    radius Ź radius+dRadius;
    FOR j: NAT IN [0..3) DO
        diff[0][j] Ź diff[0][j]+diff[1][j];
        diff[1][j] Ź diff[1][j]+diff[2][j];
        diff[2][j] Ź diff[2][j]+diff[3][j];
        ENDLOOP;
sample.value Ź sample.value+Val[];
ENDLOOP;
G3dMatrix.ReleaseMatrix[diff];
};
};