/* See SO2DPlot.h */ /* Last edited on 2003-07-08 13:33:52 by stolfi */ #include #include #include #include #include #include #include #include #include /* INTERNAL DATA TYPES */ typedef int Isoline; /* Index of an isoline. */ /* By definition, the isoline with index {k} is where the function has value {(k+0.5)*fStep}, where {fStep} is the spacing between levels. */ /* Color near zero: */ static Color ZERO_PAINT_COLOR = (Color){{1.00, 1.00, 1.00}}; /* Limiting colors for color bands (brightness = 0.5): */ static Color POS_PAINT_COLOR = (Color){{1.00, 0.33, 0.00}}; static Color NEG_PAINT_COLOR = (Color){{0.00, 0.67, 1.00}}; /* Colors for isoline drawing (brightness = 0.20): */ static Color POS_DRAW_COLOR = (Color){{0.57, 0.00, 0.00}}; static Color NEG_DRAW_COLOR = (Color){{0.00, 0.17, 1.00}}; /* Scale factor from cell-relative Y (ranges over [0_1]) to plot Y. */ //#define YSCALE (double)0.70710678118654752440 #define YSCALE (double)1.0 #define totminDepth 10 /* GLOBAL VARIABLES */ static Color *bandColor = NULL; /* color table */ /* The band between isolines {k-1} and {k} is plotted with color {bandColor[k - kPlotMin]} where {kPlotMin} is the index of the first isoline to be plotted. */ static double isolineWidth = 0.15; /* Thickness of isolines in mm. */ /* This global variable should become unnecessary once the setting of line width gets separated from the setting of line color in {ps.h}. */ /* INTERNAL PROTOTYPES */ /* The procedures {SO2DPlot_ValuesInRectangle} and {SO2DPlot_ValuesInTriangle} plot only those isolines with indices in the range {kPlotMin..kPlotMax}. The points of the root cell where {f} lies below isoline {kPlotMin} are plotted as a single color band; and ditto for values above isoline {kPlotMax}. The {layer} parameter selects between color band painting ({layer == 0}) and isoline plotting ({layer == 1}). */ void SO2DPlot_ValuesInRectangle ( SOPlot_PSStream *ps, /* Picture stream. */ SOFunction *f, /* Function to be plotted. */ FuncMap *FMap, /* Function Map for f. */ double minR[2], /* Low corner of rectangle {C} rel. root cell. */ double maxR[2], /* High corner of rectangle {C} rel. root cell. */ Isoline kPlotMin, /* Index of lowest isoline to plot. */ Isoline kPlotMax, /* Index of highest isoline to plot. */ double fStep, /* Value increment between isolines/bands. */ int layer, /* 0 = paint color bands, 1 = draw isolines. */ double *fObsMin, /* (IN/OUT) Min value seen during plot. */ double *fObsMax /* (IN/OUT) Max value seen during plot. */ ); /* Plots color bands or isolines within the rectangle with lower corner {minR} and upper corner {maxR} (in root-cell-relative coordinates). Also updates the values of {fObsMin} and {fObsMax}, which must be initialized by the caller. */ void SO2DPlot_ValuesInTriangle ( SOPlot_PSStream *ps, /* Picture stream. */ double Px, double Py, double Pf, /* Coordinates and value at vertex {P}. */ double Qx, double Qy, double Qf, /* Coordinates and value at vertex {Q}. */ double Rx, double Ry, double Rf, /* Coordinates and value at vertex {R}. */ Isoline kPlotMin, /* Index of lowest isoline to plot. */ Isoline kPlotMax, /* Index of highest isoline to plot. */ double fStep, /* Value increment between isolines/bands. */ int layer /* 0 = paint color bands, 1 = draw isolines. */ ); /* Plots color bands or isolines within the triangle with corners {(Px,Py), (Qx,Qy), (Rx,Ry)} (in root-cell-relative coordinates). Function values are interpolated between the given values {Pf,Qf,Rf} at those corners. */ void SO2DPlot_PaintRectangle ( SOPlot_PSStream *ps, /* Picture stream. */ double minR[2], double maxR[2], Color *c ); /* Paints the rectangle defined by its lower corner {minR} and its upper corner {maxR} (in root-cell-relative coordinates) with uniform color {c}. */ void SO2DPlot_PaintTriangle ( SOPlot_PSStream *ps, /* Picture stream. */ double Ax, double Ay, double Bx, double By, double Cx, double Cy, Color *c ); /* Paints the triangle defined by its corners {(Ax,Ay), (Bx,By), and (Cx,Cy)} (in root-cell-relative coordinates) with uniform color {c}. */ void SO2DPlot_DrawSegment ( SOPlot_PSStream *ps, /* Picture stream. */ double Ax, double Ay, double Bx, double By, Color *c ); /* Draws the segment with endpoints {(Ax,Ay)} and {(Bx,By)} (in root-cell-relative coordinates) with uniform color {c}. */ void SO2DPlot_SetColors ( Isoline kPlotMin, /* Index of first isoline to be drawn. */ Isoline kPlotMax, /* Index of last isoline to be drawn. */ double fStep /* Value step between isolines. */ ); /* Initializes the color table for band painting. The band between isolines {-1} and {0} is painted white. If {kPlotMin < 0}, the band below isoline {kPlotMin} is painted with {NEG_PAINT_COLOR}. If {kPlotMax > 0}, the band above isoline {kPlotMax} is painted with {POS_PAINT_COLOR}. Intermediate isolines are painted with intermediate colors. */ void SO2DPlot_Subtree ( SOPlot_PSStream *ps, /* Picture stream. */ SOGrid_Node *p, SOGrid_Index k, SOGrid_Rank rank, double minR[2], double maxR[2], int maxDepth ); /* Draws the cells of the dyadic grid rooted at node {p}, which is assumed to correspond to a cell {C} of index {k} and given {rank}. The boundary of {C} must be drawn by the client. The drawing is clipped to the rectangle with corners {minR} and {maxR} (in cell-root-relative coordinates). Cells at level higher than {maxDepth} are ignored, i.e. the subtree is implicitly truncated at level {maxDepth}. */ /* IMPLEMENTATIONS */ void SO2DPlot_AddFunctionPicture ( SOPlot_PSStream *ps, /* Picture stream. */ SOFunction *f, /* Function to be plotted. */ FuncMap *FMap, /* Function Map for f. */ SOGrid_Tree *tree, /* Reference grid for adaptive domain subdivision. */ double minR[2], /* Low corner in root-relative coordinates. */ double maxR[2], /* High corner in root-relative coordinates. */ double fPlotMin, /* Minimum value for color mapping. */ double fPlotMax, /* Maximum value for color mapping. */ double fStep, /* Value increment between isolines/bands. */ double lineWidth, /* Line width for isoline plotting. */ int minDepth, /* Minimum subdivision depth for plotting. */ int extraDepth, /* Extra depth of subdivision for leaves of {tree}. */ int maxDepth, /* Maximum subdivision depth for plotting. */ double gridWidth, /* Line width for grid drawing. */ bool grid, /* TRUE plots the tree */ int maxDepthTree, /* Omit cells below this depth for tree plotting. */ bool bands, /* TRUE paints between isolines, FALSE leaves blank. */ bool isolines, /* TRUE draws the isolines, FALSE omits them. */ double *fObsMin, /* (IN/OUT) Min value seen during plot. */ double *fObsMax, /* (IN/OUT) Max value seen during plot. */ char *caption /* Either picture comments or NULL. */ ) { double minC[2], maxC[2]; /* Coordinates of a cell {C}. */ /* Indices of first and last isolines. */ /* Allow for rounding errors in the parameters {fPlotMin}, {fmax}. */ Isoline kPlotMin = (int)(ceil(fPlotMin/fStep + 0.4999999)); Isoline kPlotMax = (int)(floor(fPlotMax/fStep + 0.5000001)); int layer; /* 0 = bands, 1 = isolines; */ /* Set isoline width for internal procs. */ isolineWidth = lineWidth; pswr_new_picture(ps, minR[0], maxR[0], minR[1], maxR[1]); auto void SO2DPlot_ValuesInCell ( SOGrid_Index k, /* The index of a dyadic (sub)cell. */ SOGrid_Rank rank, /* Its rank. */ SOGrid_Node *nd, /* A node of {tree}, or NULL. */ SOGrid_Rank lastRank /* Rank of lowest node containing the cell. */ ); /* Plots color bands (if {layer == 0} or isolines (if {layer == 1}) of function {f}, inside the dyadic cell {C} in level {rank} whose index is {k}. Assumes that {nd} is the corresponding node in {tree}, or is NULL if such node doesn't exist. In either case, assumes {lastRank} is the rank of the lowest node of {tree} whose cell contains {C}; or is 0 if {tree == NULL}. */ void SO2DPlot_ValuesInCell ( SOGrid_Index k, /* The index of a dyadic (sub)cell. */ SOGrid_Rank rank, /* Its rank. */ SOGrid_Node *nd, /* A node of {tree}, or NULL. */ SOGrid_Rank lastRank /* Rank of lowest node containing the cell. */ ) { int locMinDepth = (lastRank + extraDepth) >= totminDepth ? (lastRank + extraDepth):totminDepth; /* Get coordinates of plot cell relative to the root cell: */ SOGrid_cell_coords(k, rank, 2, minC, maxC); /* Change to absolute coordenates. */ minC[Y] *= SQRTHALF; maxC[Y] *= SQRTHALF; /* Skip if cell lies outside the desired rectangle. */ if ((minC[X] >= maxR[X]) || (maxC[X] <= minR[X])) { return; } if ((minC[Y] >= maxR[Y]) || (maxC[Y] <= minR[Y])) { return; } /* Subdivide and recurse, or plot undivided: */ if (((rank < maxDepth)||(rank < totminDepth)) && ((rank < minDepth)||(rank < locMinDepth))) { /* Split {C} into subcells, even if it is a leaf cell. */ SOGrid_Node *nd0 = (nd ? nd->c[0]: NULL); SOGrid_Node *nd1 = (nd ? nd->c[1]: NULL); SO2DPlot_ValuesInCell(2*k, rank+1, nd0, lastRank + (nd0 != NULL)); SO2DPlot_ValuesInCell(2*k+1, rank+1, nd1, lastRank + (nd1 != NULL)); } else { /* Clip cell to given rectangle: */ if (maxC[X] > maxR[X]) { maxC[X] = maxR[X]; } if (minC[X] < minR[X]) { minC[X] = minR[X]; } if (maxC[Y] > maxR[Y]) { maxC[Y] = maxR[Y]; } if (minC[Y] < minR[Y]) { minC[Y] = minR[Y]; } /* Plot approximate isolines in rectangle. */ SO2DPlot_ValuesInRectangle ( ps, f, FMap, minC, maxC, kPlotMin, kPlotMax, fStep, layer, fObsMin, fObsMax ); } // DEBUG fprintf (stderr, "-InCell cell = %ld rank = %d\n", k, rank); } // DEBUG fprintf (stderr, " minDepth = %d extraDepth = %d maxDepth = %d\n", minDepth, extraDepth, maxDepth); if (! (isolines || bands)) { return; } assert(fPlotMin < fPlotMax, "bad function range"); assert(kPlotMin <= kPlotMax, "no isolines in range"); SO2DPlot_SetColors(kPlotMin, kPlotMax, fStep); if (bands) { layer = 0; SO2DPlot_ValuesInCell (1, 0, (tree ? tree->root : NULL), 0); } if (isolines) { layer = 1; SO2DPlot_ValuesInCell (1, 0, (tree ? tree->root : NULL), 0); } fprintf(stderr, "SO2DPlot_Values: fObsMin = %.16g, fObsMax = %.16g\n", (*fObsMin), (*fObsMax) ); if (grid && (tree != NULL)) { SO2DPlot_Tree(ps, tree, minR, maxR, gridWidth, maxDepthTree); } pswr_set_pen(ps, 0.0, 0.0, 0.0, lineWidth, 0.0, 0.0); if(caption != NULL)pswr_add_caption(ps, caption, 0.0); } /* Number of sample points in each rectangle (either 4 or 8): */ #define NSAMPLES 8 void SO2DPlot_ValuesInRectangle ( SOPlot_PSStream *ps, /* Picture stream. */ SOFunction *f, /* Function to be plotted. */ FuncMap *FMap, /* Function Map for f. */ double minR[2], /* Low corner of rectangle {C} rel. root cell. */ double maxR[2], /* High corner of rectangle {C} rel. root cell. */ Isoline kPlotMin, /* Min isoline index to plot. */ Isoline kPlotMax, /* Max isoline index to plot. */ double fStep, /* Value increment between isolines/bands. */ int layer, /* 0 = paint color bands, 1 = draw isolines. */ double *fObsMin, /* (IN/OUT) Min value seen during plot. */ double *fObsMax /* (IN/OUT) Max value seen during plot. */ ) { double p[(NSAMPLES+1)*2]; /* Sampling points, root-relative coords. */ /* Samples {[0..NSAMPLES-1]} lie around the perimeter of {C}, */ /* in CCW order; sample {[NSAMPLES]} is the center of the rectangle. */ double fp[NSAMPLES+1]; /* Sampled function values: {fp[i] = f(p[i])}. */ double *ctr = &(p[NSAMPLES*2]); double fpMin, fpMax; /* Min and max of {fp[0..NSAMPLES(either 4 or 8)]}. */ Isoline kpMin, kpMax; /* Isolines in {[fpMin _ fpMax]} are {kpMin..kmax}. */ int i; // DEBUG fprintf (stderr, "+InRectangle\n"); /* Coordinates of rectangle center: */ ctr[X] = (minR[X]+maxR[X])/2.0; ctr[Y] = (minR[Y]+maxR[Y])/2.0; /* Gets coordinates and function values of corners and mid-edges: */ /* Important: {fp[0..7]} must be in counterclockwise order. */ for (i = 0; i < NSAMPLES; i++) { double *pi = &(p[i*2]); if (NSAMPLES == 4) { /* Samples are corners only. */ pi[X] = ((i==0) || (i==3) ? minR[X] : maxR[X]); pi[Y] = ((i==0) || (i==1) ? minR[Y] : maxR[Y]); } else { /* Samples are corners and mid-edges. */ if (i % 2 == 0) { /* Corner number {i/2}: */ pi[X] = ((i==0) || (i==6) ? minR[X] : maxR[X]); pi[Y] = ((i==0) || (i==2) ? minR[Y] : maxR[Y]); } else { /* Mid-edge between corners {i} and {i+1}: */ pi[X] = ((i==1) || (i==5) ? ctr[X] : (i==7 ? minR[X] : maxR[X])); pi[Y] = ((i==7) || (i==3) ? ctr[Y] : (i==1 ? minR[Y] : maxR[Y])); } } /* Pull sample point slightly inside the cell, for discntinuous splines: */ pi[X] = 0.9999999 * pi[X] + 0.0000001 * ctr[X]; pi[Y] = 0.9999999 * pi[Y] + 0.0000001 * ctr[Y]; /* Evaluate function at sample point: */ f->m->eval(f, pi, &fp[i]); if(FMap->map != NULL) { double fmv[FMap->vDim]; int mapi; FMap->map(&fp[i], pi, fmv); for(mapi = 0; mapi < FMap->uDim; mapi++) fp[i+mapi]=fmv[mapi]; } // printf(" ** SO2DPlot DEBUG fp[i]: %.16g, i: %d \n\n", fp[i], i); // DEBUG fprintf (stderr, " f(%8.6f,%8.6f) = %10.6f\n", pi[X], pi[Y], fp[i]); } /* The center value: */ f->m->eval(f, ctr, &fp[NSAMPLES]); if(FMap->map != NULL) { double fmv[FMap->vDim]; int mapi; FMap->map(&fp[NSAMPLES], ctr, fmv); for(mapi = 0; mapi < FMap->uDim; mapi++) fp[NSAMPLES+mapi]=fmv[mapi]; } // printf(" ** SO2DPlot DEBUG fp[NSAMPLES]: %.16g, NSAMPLES: %d \n\n", fp[NSAMPLES], NSAMPLES); // DEBUG fprintf (stderr, " f(%8.6f,%8.6f) = %10.6f\n", ctr[X], ctr[Y], fp[NSAMPLES]); /* Find min and max among those sampled values: */ fpMin = fp[0]; fpMax = fp[0]; for(i = 1; i < NSAMPLES+1; i++) { if (fp[i] < fpMin) { fpMin = fp[i]; } if (fp[i] > fpMax) { fpMax = fp[i]; } } /* Update global min and max: */ if (fpMin < (*fObsMin)) { (*fObsMin) = fpMin; } if (fpMax > (*fObsMax)) { (*fObsMax) = fpMax; } /* Compute first and last isolines that enter this rectangle: */ kpMin = (int)(ceil(fpMin/fStep - 0.4999999)); kpMax = (int)(floor(fpMax/fStep - 0.5000001)); /* Do any isolines enter the rectangle? */ if ((kpMin > kPlotMax) || (kpMax < kPlotMin) || (kpMin > kpMax)) { /* No. */ if (layer == 0) { /* Paint whole rectangle with uniform color. */ Color *cBand; int kBand = kpMin; if (kpMin > kPlotMax) { kBand = kPlotMax + 1; } if (kpMax < kPlotMin) { kBand = kPlotMin; } cBand = &(bandColor[kBand - kPlotMin]); SO2DPlot_PaintRectangle(ps, minR, maxR, cBand); } } else { /* Paint each triangle: */ for (i = 0; i < NSAMPLES; i++) { int j = (i + 1) % NSAMPLES; double *pi = &(p[2*i]); double *pj = &(p[2*j]); SO2DPlot_ValuesInTriangle ( ps, ctr[X], ctr[Y], fp[NSAMPLES], pi[X], pi[Y], fp[i], pj[X], pj[Y], fp[j], kPlotMin, kPlotMax, fStep, layer ); } } // DEBUG fprintf (stderr, "-InRectangle\n"); } void SO2DPlot_ValuesInTriangle ( SOPlot_PSStream *ps, /* Picture stream. */ double Px, double Py, double Pf, double Qx, double Qy, double Qf, double Rx, double Ry, double Rf, Isoline kPlotMin, Isoline kPlotMax, double fStep, int layer ) { double temp; Isoline kMin, kMax; /* Isolines that enter the triangle are {kMin..kMax}. */ // DEBUG fprintf (stderr, "+InTriangle\n"); /* Permute corners so that {Pf <= Qf <= Rf}: */ if( Pf > Qf ) { temp = Qx; Qx = Px; Px = temp; temp = Qy; Qy = Py; Py = temp; temp = Qf; Qf = Pf; Pf = temp; } if( Qf > Rf ) { temp = Qx; Qx = Rx; Rx = temp; temp = Qy; Qy = Ry; Ry = temp; temp = Qf; Qf = Rf; Rf = temp; } if( Pf > Qf ) { temp = Qx; Qx = Px; Px = temp; temp = Qy; Qy = Py; Py = temp; temp = Qf; Qf = Pf; Pf = temp; } assert( Pf <= Qf, "Pf > Qf" ); assert( Qf <= Rf, "Qf > Rf" ); /* Compute first and last isolines that enter this rectangle: */ kMin = (int)(ceil(Pf/fStep - 0.4999999)); kMax = (int)(floor(Rf/fStep - 0.5000001)); // DEBUG fprintf(stderr, " Pf = %10.8f kMin = %d Rf = %10.8f kMax = %d\n", Pf, kMin, Rf, kMax); /* Do any isolines enter the triangle? */ if ((kMin > kPlotMax) || (kMax < kPlotMin) || (kMin > kMax)) { /* No. */ // DEBUG fprintf(stderr, " kPlotMin = %d kPlotMax = %d\n", kPlotMin, kPlotMax); if (layer == 0) { /* Paint whole triangle with color: */ Color *cBand; int kBand = kMin; if (kMin > kPlotMax) { kBand = kPlotMax + 1; } if (kMax < kPlotMin) { kBand = kPlotMin; } cBand = &(bandColor[kBand - kPlotMin]); SO2DPlot_PaintTriangle(ps, Px, Py, Qx, Qy, Rx, Ry, cBand); } } else { /* Break the triangle {P,Q,R} into slices at isoline levels. */ double Lf = Pf - 1.0; /* Function value of previous isoline. */ double LUx = Px, LUy = Py; /* Lower point on seg {P-R}. */ double LVx = Px, LVy = Py; /* Lower point on seg {P-Q} or {Q-R}. */ bool Lsame = TRUE; /* Do {LU} and {LV} coincide? */ Isoline k; /* Restrict to specified isoline range: */ if (kMin < kPlotMin) { kMin = kPlotMin; } if (kMax > kPlotMax) { kMax = kPlotMax; } for (k = kMin; k <= kMax + 1; k++) { /* Function value at higher isoline: */ double Hf = (k > kMax ? Rf + 1.0 : (k + 0.5)*fStep); double HUx, HUy; /* Higher point on seg {P-R}. */ double HVx, HVy; /* Higher point on seg {P-Q} or {Q-R}. */ bool Hsame; /* Do {HU} and {HV} coincide? */ /* Slice index {k} is defined by: two lower points {LU = (LUx,LUy)} and {LV = (LVx,LVy)}; two higher points {HU = (HUx,HUy)} and {HV = (HVx,HVy)}; and possibly the point {Q}, if {Lf= Pf, "Hf < Pf"); // DEBUG fprintf(stderr, " k = %d Hf = %10.6f\n", k, Hf); /* Compute points {HU,HV} on upper isoline: */ if (Hf >= Rf) { /* Last band: */ HUx = HVx = Rx; HUy = HVy = Ry; Hsame = TRUE; } else { double sU = (Hf - Pf)/(Rf - Pf), tU = 1.0 - sU; HUx = tU * Px + sU * Rx; HUy = tU * Py + sU * Ry; if (Hf <= Qf) { double sV = (Hf - Pf)/(Qf - Pf), tV = 1.0 - sV; HVx = tV * Px + sV * Qx; HVy = tV * Py + sV * Qy; } else { double sV = (Hf - Qf)/(Rf - Qf), tV = 1.0 - sV; HVx = tV * Qx + sV * Rx; HVy = tV * Qy + sV * Ry; } Hsame = FALSE; } if (layer == 1) { /* Draw the isoline at level {k-1}, if it is non-empty. */ if ((k > kMin) && (Lf > Pf)) { Color *c = (k > 0 ? &(POS_DRAW_COLOR) : &(NEG_DRAW_COLOR)); SO2DPlot_DrawSegment(ps, LUx, LUy, LVx, LVy, c); } } else { /* Split the slice into trianglets, and paint them: */ Color *c = &(bandColor[k - kPlotMin]); if ((Lf >= Qf) || (Hf <= Qf)) { /* Region is the trapezoid {LU,LV,HV,HU}. */ if (Lsame && Hsame) { /* Nothing to do */ } else if (Lsame) { SO2DPlot_PaintTriangle(ps, LUx, LUy, HVx, HVy, HUx, HUy, c); } else if (Hsame) { SO2DPlot_PaintTriangle(ps, LUx, LUy, LVx, LVy, HVx, HVy, c); } else { SO2DPlot_PaintTriangle(ps, LUx, LUy, HVx, HVy, HUx, HUy, c); SO2DPlot_PaintTriangle(ps, HVx, HVy, LUx, LUy, LVx, LVy, c); } } else { /* Region is the convex pentagon {LU,LV,Q,HV,HU}. */ if (Lsame && Hsame) { SO2DPlot_PaintTriangle(ps, LUx, LUy, Qx, Qy, HUx, HUy, c); } else if (Lsame) { SO2DPlot_PaintTriangle(ps, LUx, LUy, Qx, Qy, HVx, HVy, c); SO2DPlot_PaintTriangle(ps, LUx, LUy, HVx, HVy, HUx, HUy, c); } else if (Hsame) { SO2DPlot_PaintTriangle(ps, LUx, LUy, LVx, LVy, HVx, HVy, c); SO2DPlot_PaintTriangle(ps, LVx, LVy, Qx, Qy, HVx, HVy, c); } else { SO2DPlot_PaintTriangle(ps, LVx, LVy, Qx, Qy, LUx, LUy, c); SO2DPlot_PaintTriangle(ps, LUx, LUy, Qx, Qy, HUx, HUy, c); SO2DPlot_PaintTriangle(ps, HUx, HUy, Qx, Qy, HVx, HVy, c); } } } LUx = HUx; LUy = HUy; LVx = HVx; LVy = HVy; Lf = Hf; Lsame = Hsame; } } // pswr_set_pen(ps, 0.0, 0.0, 0.0, 0.10, 0.0,0.0); // pswr_segment(ps, Px, Py*YSCALE, Qx, Qy*YSCALE); // pswr_segment(ps, Qx, Qy*YSCALE, Rx, Ry*YSCALE); // pswr_segment(ps, Rx, Ry*YSCALE, Px, Py*YSCALE); // DEBUG fprintf (stderr, "-InTriangle\n"); } void SO2DPlot_PaintRectangle ( SOPlot_PSStream *ps, /* Picture stream. */ double minR[2], double maxR[2], Color *c ) { double R = c->c[0], G = c->c[1], B = c->c[2]; if ((R >= 0) && (G >= 0) && (B >= 0)) { double Lx = minR[X], Ly = minR[Y] * YSCALE; double Hx = maxR[X], Hy = maxR[Y] * YSCALE; pswr_set_fill_color(ps, R, G, B); pswr_rectangle(ps, Lx, Hx, Ly, Hy, 1, 0); } } void SO2DPlot_PaintTriangle ( SOPlot_PSStream *ps, /* Picture stream. */ double Ax, double Ay, double Bx, double By, double Cx, double Cy, Color *c ) { double R = c->c[0], G = c->c[1], B = c->c[2]; if ((R >= 0) && (G >= 0) && (B >= 0)) { Ay *= YSCALE; By *= YSCALE; Cy *= YSCALE; pswr_set_fill_color(ps, R, G, B); pswr_triangle(ps, Ax, Ay, Bx, By, Cx, Cy, 1, 0); /* pswr_set_pen(ps, 0.5, 0.5, 0.7, 0.001, 0.0,0.0); */ /* pswr_segment(ps, Ax, Ay, Bx, By); */ /* pswr_segment(ps, Ax, Ay, Cx, Cy); */ /* pswr_segment(ps, Bx, By, Cx, Cy); */ } } void SO2DPlot_DrawSegment ( SOPlot_PSStream *ps, /* Picture stream. */ double Ax, double Ay, double Bx, double By, Color *c ) { double R = c->c[0], G = c->c[1], B = c->c[2]; if ((R >= 0) && (G >= 0) && (B >= 0)) { pswr_set_pen(ps, R, G, B, isolineWidth, 0.0,0.0); pswr_segment(ps, Ax, Ay*YSCALE, Bx, By*YSCALE); } } void SO2DPlot_SetColors ( Isoline kPlotMin, /* Index of first isoline to be drawn. */ Isoline kPlotMax, /* Index of last isoline to be drawn. */ double fStep /* Value step between isolines. */ ) { int nBands = (kPlotMax + 1) - kPlotMin + 1; int i; if (bandColor != NULL) { free(bandColor); } bandColor = (Color *)malloc(nBands*sizeof(Color)); assert(bandColor != NULL, "No memory for color table allocation"); for(i = 0; i < nBands; i++) { Isoline k = kPlotMin + i; /* Index of isoline above band. */ if (k < 0) { /* Values are negative. */ double s = ((double)k)/((double)kPlotMin); bandColor[i] = SOPlot_ColorScale(s, &ZERO_PAINT_COLOR, &NEG_PAINT_COLOR); } else if (k > 0) { /* Values are positive. */ double s = ((double)k)/((double)kPlotMax); bandColor[i] = SOPlot_ColorScale(s, &ZERO_PAINT_COLOR, &POS_PAINT_COLOR); } else { bandColor[i] = ZERO_PAINT_COLOR; } } } void SO2DPlot_Tree( SOPlot_PSStream *ps, /* Picture stream. */ SOGrid_Tree *tree, /* The grid to plot. */ double minR[2], /* Low corner of region to plot. */ double maxR[2], /* High corner of region to plot. */ double lineWidth, /* Line width for grid cell boundaries. */ int maxDepth /* Omit cells below this depth. */ ) { /* pswr_set_pen(ps, 0.0, 0.0, 0.0, lineWidth, 0.0, 0.0); */ pswr_set_pen(ps, 0.6, 0.6, 0.6, lineWidth, 0.0, 0.0); pswr_rectangle(ps, minR[X], maxR[X], YSCALE*minR[Y], YSCALE*maxR[Y], 0, 1); if (tree == NULL) { return; } assert(tree->d == 2, "wrong tree dimension"); SO2DPlot_Subtree(ps, tree->root, 1, 0, minR, maxR, maxDepth); } void SO2DPlot_Subtree ( SOPlot_PSStream *ps, /* Picture stream. */ SOGrid_Node *p, SOGrid_Index k, SOGrid_Rank rank, double minR[2], double maxR[2], int maxDepth ) { double minC[2], maxC[2]; /* Ignore cells at levels below {maxDepth}: */ if (rank >= maxDepth) { return; } /* If cell is childless, there is nothing to do: */ if ((p == NULL) || ((p->c[0] == NULL) && (p->c[1] == NULL))) { return; } SOGrid_cell_coords(k, rank, 2, minC, maxC); /* Change to absolute coordenates. */ minC[Y] *= SQRTHALF; maxC[Y] *= SQRTHALF; /* Skip plot if cell lies outside the region of interest: */ if ((minC[X] >= maxR[X]) || (maxC[X] <= minR[X])) { return; } if ((minC[Y] >= maxR[Y]) || (maxC[Y] <= minR[Y])) { return; } /* Draw the interior of daughter cells first: */ SO2DPlot_Subtree(ps, p->c[0], 2*k, rank + 1, minR, maxR, maxDepth); SO2DPlot_Subtree(ps, p->c[1], 2*k+1, rank + 1, minR, maxR, maxDepth); /* Now draw the boundary between the two daughter cells: */ if ((rank % 2) == 0) { /* Longest axis is {X}: */ double xMid = (minC[X]+maxC[X])/2; double yMin = YSCALE * (minC[Y] < minR[Y] ? minR[Y] : minC[Y]); double yMax = YSCALE * (maxC[Y] > maxR[Y] ? maxR[Y] : maxC[Y]); pswr_segment(ps, xMid, yMin, xMid, yMax); } else { /* Longest axis is {Y}: */ double yMid = YSCALE * (minC[Y]+maxC[Y])/2; double xMin = (minC[X] < minR[X] ? minR[X] : minC[X]); double xMax = (maxC[X] > maxR[X] ? maxR[X] : maxC[X]); pswr_segment(ps, xMin, yMid, xMax, yMid); } } void SO2DPlot_SingleFunction ( bool eps, /* TRUE for EPS figures, FALSE for PS document. */ char *paperSize, /* Paper size ({letter}, {a3}, etc.). */ double hFigSize, /* Total figure width (in pt). */ double vFigSize, /* Total figure height (in pt). */ SOFunction *f, /* Function to plot. */ FuncMap *FMap, /* Function Map for f. */ SOGrid_Tree *tree, /* Optional grid to plot. */ double minR[2], /* Low corner in root-relatives coordinates. */ double maxR[2], /* High corner in root-relative coordinates. */ double fPlotMin, /* Nominal minimum {f} value, for color scale. */ double fPlotMax, /* Nominal maximum {f} value, for color scale. */ double fStep, /* Value step for isolines and color bands. */ double isolineWidth, /* Line width for isoline plotting. */ double gridWidth, /* Line width for grid drawing. */ int meshDepth, /* Depth of bisection recursion for plotting. */ bool bands, /* TRUE plots color bands. */ bool isolines, /* TRUE draws isolines. */ bool grid, /* TRUE plots the tree */ char *figName, /* Figure name (minus extension). */ double *fObsMin, /* (IN/OUT) Min value seen during plot. */ double *fObsMax /* (IN/OUT) Max value seen during plot. */ ) { int maxDepthIsoline = meshDepth; int minDepthIsoline = meshDepth - 6; int extraDepthIsoline = 2; int maxDepthTree = meshDepth + 20; double hPageSize, vPageSize; /* printf("DEBUG - hFigSize: %f, vFigSize: %f, paperSize: %s \n", hFigSize, vFigSize, paperSize); */ if(paperSize == NULL) {hPageSize = hFigSize + 150; vPageSize = (vFigSize + 150)*SQRTHALF;} else pswr_get_paper_dimensions(paperSize, &hPageSize, &vPageSize); /* printf("DEBUG - hPageSize: %f, vPageSize: %f, paperSize: %s \n", hPageSize, vPageSize, paperSize); */ SOPlot_PSStream *ps = SOPlot_OpenStream(eps, figName, paperSize, hPageSize, vPageSize); double Margin = (eps ? 5 : 50); pswr_set_page_layout(ps, hPageSize - (2*Margin), vPageSize - (2*Margin), 0, Margin, Margin, (! eps), 1, 1); if (isolines || bands) { SO2DPlot_AddFunctionPicture ( ps, f, FMap, tree, minR, maxR, fPlotMin, fPlotMax, fStep, isolineWidth, minDepthIsoline, extraDepthIsoline, maxDepthIsoline, gridWidth, grid, maxDepthTree, bands, isolines, fObsMin, fObsMax, figName ); } SOPlot_CloseStream(ps); } /* Limits to avoid excessive plotting for bad choices of {fStep}: */ /* Should make them consistent with internal limits of {SO2DPlot_Values}. */ #define DefaultIsolines 10 #define MaxIsolinesInRange (3*DefaultIsolines/2) void SO2DPlot_FixRange ( double *fRange, double *fStep, double alpha, SOFunction *f, int verbose ) { /* If no place to return result, do nothing: */ if (fRange == NULL) { return; } mumble(1, "checking nominal function range...\n"); /* If given range is invalid, ignore it: */ if ((*fRange) <= 0.0) { mumble(1, " given range is invalid, ignoring it.\n"); alpha = 1.0; } /* If no adjustment allowed, do nothing: */ if (alpha != 0.0) { double newRange; { /* Estimate the maximum function value: */ double estRange = SO2DPlot_EstRange(f, verbose); if (alpha == 1.0) { newRange = estRange; } else { double oldRange = *fRange; mumble(1, " mixing with given range (alpha = %8.6f)...\n", alpha); newRange = alpha * estRange + (1.0 - alpha)*oldRange; } } assert(newRange > 0.0, "empty range - something failed"); mumble(1, " rounding computed range (%10.4e) to a nice value...\n", newRange); { double step = SO2DPlot_RoundToNice(newRange/DefaultIsolines); newRange = DefaultIsolines * step; *fRange = newRange; mumble(1, " new range is [ %10.4e _ %10.4e ].\n", -(*fRange), (*fRange)); } } /* Fix {*fStep} if needed: */ if (fStep == NULL) { return; } mumble(1, "checking isoline spacing...\n"); if ((*fStep) <= (*fRange)/MaxIsolinesInRange) { (*fStep) = SO2DPlot_RoundToNice((*fRange)/DefaultIsolines); mumble(1, " given spacing is too small, fixed to %10.4e.\n", *fStep); } } double SO2DPlot_RoundToNice(double x) { double z = 1.0; /* Allow for rounding errors in the algorithm below: */ x = 0.9999999 * x; /* If is too small, don't bother: */ if (x < 1.0e-300) { return 1.0e-300; } /* Finds the smallest power of 10 that is no less than {x}: */ while (z/10.0 > x) { z /= 10.0; } while (z <= x ) { z *= 10.0; } /* Finds a nice multiple of that power of 10: */ if (x < 0.25 * z) { return 0.25 * z; } else if (x < 0.5 * z) { return 0.5 * z; } else { return z; } } double SO2DPlot_EstRange(SOFunction *f, int verbose) { assert(FALSE, "not implemented yet"); return 0.0; // SOSpline *fpw; // // auto double func(DomPoint *p); // // double func(DomPoint *p) // { return f->m->eval(f, p); } // // if ((fpw = SOSpline_Cast(f)) != NULL) // { mumble(1, " estimating max function value (grid sampling)...\n"); // return SORange_OnSOGrid(func, fpw->d->tree, 40); // } // else // { mumble(1, " estimating max function value (root cell sampling)...\n"); // return SORange_OnRootCell(func, 40) // } }