/* Last edited on 2024-12-29 00:00:19 by stolfi */ /* ---------------------------------------------------------------- */ /* Replaced by Rafael's latest code: */ /* The parameter {avgWidth} specifies the width of the kernel used when reducing the slope map and expanding the height map. See {pst_slope_map_shrink} and {pst_height_map_expand}. The weight map is reduced by arithmetic averagin if {harmAvgW} is false, by harmonic averaging if {harmAvgW} is true (see {pst_weight_map_shrink}). If {newStyle} is true the recursion stops at 1x1 pixel */ float_image_t *pst_height_map_expand(float_image_t *JZ, int NXI, int NYI, int expOrder) { /* !!! Generalize and move to {float_image_mscale.h} !!! */ demand(expOrder == 2, "{expOrder} other than 2 is not supported yet"); assert(JZ->sz[0] == 1); int NXJ = JZ->sz[1]; assert(NXI/2 + 1 == NXJ); int NYJ = JZ->sz[2]; assert(NYI/2 + 1 == NYJ); float_image_t *IZ = float_image_new(1, NXI, NYI); int xI, yI; for(yI = 0; yI < NYI; yI++) { for(xI = 0; xI < NXI; xI++) { /* Get relevant samples from original image: */ int xJ0 = xI/2, xJ1 = (xI+1)/2; int yJ0 = yI/2, yJ1 = (yI+1)/2; assert(xJ1 < NXJ); assert(yJ1 < NYJ); double d00 = float_image_get_sample(JZ, 0, xJ0, yJ0); double d01 = (yJ1 == yJ0 ? d00 : float_image_get_sample(JZ, 0, xJ0, yJ1)); double d10 = (xJ1 == xJ0 ? d00 : float_image_get_sample(JZ, 0, xJ1, yJ0)); double d11 = ((xJ1 == xJ0) && (yJ1 == yJ0) ? d00 : float_image_get_sample(JZ, 0, xJ1, yJ1)); /* Average samples, and scale {Z} up by 2: */ float v = (d00 + d01 + d10 + d11)/2; /* Store sample in expanded image: */ float_image_set_sample(IZ, 0, xI, yI, v); } } return IZ; } float_image_t *pst_slope_map_shrink_old(float_image_t *IG, float_image_t *IW); /* Older version of {pst_slope_map_shrink}. !!! Delete after the new one is tested. !!! */ float_image_t *pst_slope_map_shrink(float_image_t *IG, float_image_t *IW, int avgWidth); float_image_t *pst_slope_map_shrink(float_image_t *IG, float_image_t *IW, int avgWidth) { int NX_JG = (IG->sz[1]+1)/2; int NY_JG = (IG->sz[2]+1)/2; int dxy = (avgWidth-1)/2; return float_image_mscale_shrink(IG, IW, NX_JG, NY_JG, dxy, dxy, avgWidth); } float_image_t *pst_slope_map_shrink_old(float_image_t *IG, float_image_t *IW) { int NC = IG->sz[0]; int NXI = IG->sz[1]; int NXJ = (NXI+1)/2; int NYI = IG->sz[2]; int NYJ = (NYI+1)/2; if (IW != NULL) { assert(IW->sz[0] == 1); assert(IW->sz[1] == NXI); assert(IW->sz[2] == NYI); } float_image_t *JG = float_image_new(NC, NXJ, NYJ); int xJ, yJ; for(yJ = 0; yJ < NYJ; yJ++) { for(xJ = 0; xJ < NXJ; xJ++) { int xI = 2*xJ, yI = 2*yJ; double w00 = 1.0, w01 = 1.0, w10 = 1.0, w11 = 1.0; /* Weights of each sample: */ if (IW == NULL) { w00 = get_sample(IW, 0, xI, yI); w10 = ((xI+1 < NXI) ? get_sample(IW, 0, xI+1, yI) : w00); w01 = ((yI+1 < NYI) ? get_sample(IW, 0, xI, yI+1) : w00); w11 = ((xI+1 < NXI) && (yI+1 < NYI) ? get_sample(IW, 0, xI+1, yI+1) : w00); } double wtot = w00 + w01 + w10 + w11; if (wtot == 0.0) { w00 = w01 = w10 = w11 = 0.25; wtot = 1.0; } int c; for (c = 0; c < NC; c++) { /* Samples in original image (with border padding): */ double d00 = get_sample(IG, c, xI, yI); double d10 = ((xI+1 < NXI) ? get_sample(IG, c, xI+1, yI) : d00); double d01 = ((yI+1 < NYI) ? get_sample(IG, c, xI, yI+1) : d00); double d11 = ((xI+1 < NXI) && (yI+1 < NYI) ? get_sample(IG, c, xI+1, yI+1) : d00); float v = (w00*d00 + w01*d01 + w10*d10 + w11*d11)/wtot; /* All three coordinate {X,Y,Z} get scaled, so the derivatives don't: */ /* Set sample in shrunken image: */ set_sample(JG, c, xJ, yJ, v); } } } return JG; } imgsys_t* pst_slope_map_build_integration_system(float_image_t *G, float_image_t *W, bool_t full) { /* Get/check the sizes of the slope maps: */ int NC = G->sz[0]; assert(NC == 2); int NX_G = G->sz[1]; int NY_G = G->sz[2]; if (W != NULL) { assert(W->sz[0] == 1); assert(W->sz[1] == NX_G); assert(W->sz[2] == NY_G); } /* Check the size of the system: */ int NX_Z = NX_G+1; int NY_Z = NY_G+1; int NXY_Z = NX_Z*NY_Z; /* Number of unknowns (heights). */ /* Conceptually, we define two quadratic mismatch terms {qx[x,y](h)} and {qy[x,y](h)} for each edge in the pixel grid. The term {qx[x,y](h)} is the square of the difference between the mean slope along the edge from {(x,y)} to {(x+1,y)}, as estimated from the slope map {G}, and the mean slope computed by numerical difference of the (unknown) heights at the two endpoints of that edge. This terms exists for all {x} in {0..NX_G-1} and all {y} in {0..NY_G}. { qx[x,y](h) = (dx[x,y] - (- h[k00] + h[k10]))^2 } where {h[k00]} and {h[k10]} are the {h} variables corresponding to heights {z[x,y]} and {z[x+1,y]}, respectively; and {dx} is the interpolated slope { dx[x,y] = (wpm*dxpm + wpp*dxpp)/wx[x,y] } where {wpm == IW[0,x,y-1]}, {wpp == IW[0,x,y]}, {wx[x,y] == wpm+wpp}, {dxpm == G[0,x,y-1]}, and {dxpp == G[0,x,y]}. The term {qy[x,y](h)} is similarly defined for edges parallel to the Y axis. It exists for all {x} in {0..NX_G} and all {y} in {0..NY_G-1}. { qy[x,y](h) = (dy[x,y] - (- h[k00] + h[k01]))^2 } { dy[x,y] = (wmp*dxmp + wpp*dxpp)/wy[x,y] } where {wmp == IW[0,x,y-1]}, {wy[x,y] == wmp+wpp}, {dymp == G[1,x-1,y]}, and {dypp == G[1,x,y]}. We want to minimize the sum {Q(h)} of all terms {qx[x,y](h)} and {qy[x,y](h)}, each weighted by the corresponding edge weight {wx[x,y]} or {wy[x,y]}. Note that if the slope map {G} is zero, we can anihilate all those terms (and therefore Q) by setting all heights equal to 1. Therefore, if the system is to be solved iteratively, one must take care to exclude this solution. Minimizing {Q} is equivalent to finding the {h} vector that anihilates the gradient {G} of {Q}. The gradient of {Q} is an affine function of {h}, namely {G(h) = A h - b}. Each term {wx[x,y]*qx[x,y](h)} contributes to two components of the gradient, namely those with indices {k00,k10}. Similarly, each term {wy[x,y]*qy[x,y]} contributes to two components, with indices {k00,k01}. Thus, to build the system, we add those contributions to the appropriate elements of {A} and {b}. */ /* Make sure that we can build the system: */ assert(MAXCOEFS >= 5); /* Count the elements that are non-null and assemble the {h}-index table {ix[0..NXY_Z-1]}: */ long int *ix = (long int *)notnull(malloc(NXY_Z*sizeof(long int)), "no mem"); long int N; if (full || (W == NULL)) { N = NXY_Z; long int xy; for (xy = 0; xy < NXY_Z; xy++) { ix[xy] = xy; } } else { N = 0; int x, y; for(y = 0; y <= NY_G; y++) { for(x = 0; x <= NX_G; x++) { /* Compute the position in the {ix} array: */ long int xy = x + y*NX_Z; /* Get the slopes and relative weights of the pixels incident to {(x,y)}: */ bool_t determined = ((x > 0) && (y > 0) && get_sample(W, 0, x-1, y-1) > 0) || ((x > 0) && (y < NY_G) && get_sample(W, 0, x-1, y) > 0) || ((x < NX_G) && (y > 0) && get_sample(W, 0, x, y-1) > 0) || ((x < NX_G) && (y < NY_G) && get_sample(W, 0, x, y) > 0); if (determined) { /* Some neighboring slope has positive weight. */ /* Add the element {Z[x,y]} to the system: */ ix[xy] = N; N++; } else { /* All neighboring slopes have weight 0. */ /* Exclude the element {Z[x,y]} from the system: */ ix[xy] = -1; } } } } assert(N <= NXY_Z); /* Build the inverse tables {col[0..N-1],row[0..N-1]}: */ int *col = (int *)notnull(malloc(N*sizeof(int)), "no mem"); int *row = (int *)notnull(malloc(N*sizeof(int)), "no mem"); { long int xy; for (xy = 0; xy < NXY_Z; xy++) { long int k = ix[xy]; if (k >= 0) { col[k] = xy%NX_Z; row[k] = xy/NX_Z; } } } /* Allocate the system's equations: */ imgsys_t *S = pst_imgsys_new(NX_Z, NY_Z, N, ix, col, row); auto long int ind_h(int x, int y); /* This function computes the index into the solution vector {h} corresponding to vertex {[x,y]} of the images. Returns -1 if {Z[x,y]} is not in the system. */ long int ind_h(int x, int y) { if ((x < 0) || (x >= NX_Z)) { return -1; } if ((y < 0) || (y >= NY_Z)) { return -1; } return S->ix[x + y*NX_Z]; } /* Now build the equations for the variables in the system: */ long int k; for(k = 0; k < N; k++) { /* Get the indices of the pixel corresponding to the unknown {h[k]}: */ int x = col[k]; int y = row[k]; /* Get the slopes and relative weights of the pixels incident to {(x,y)}: */ double dxmm = 0.0, dxmp = 0.0, dxpm = 0.0, dxpp = 0.0; double dymm = 0.0, dymp = 0.0, dypm = 0.0, dypp = 0.0; double wmm = 0.0, wmp = 0.0, wpm = 0.0, wpp = 0.0; if ((x > 0) && (y > 0)) { dxmm = get_sample(G, 0, x-1, y-1); dymm = get_sample(G, 1, x-1, y-1); wmm = (W == NULL ? 1.0 : get_sample(W, 0, x-1, y-1)); assert(wmm >= 0.0); } if ((x > 0) && (y < NY_G)) { dxmp = get_sample(G, 0, x-1, y); dymp = get_sample(G, 1, x-1, y); wmp = (W == NULL ? 1.0 : get_sample(W, 0, x-1, y)); assert(wmp >= 0.0); } if ((x < NX_G) && (y > 0)) { dxpm = get_sample(G, 0, x, y-1); dypm = get_sample(G, 1, x, y-1); wpm = (W == NULL ? 1.0 : get_sample(W, 0, x, y-1)); assert(wpm >= 0.0); } if ((x < NX_G) && (y < NY_G)) { dxpp = get_sample(G, 0, x, y); dypp = get_sample(G, 1, x, y); wpp = (W == NULL ? 1.0 : get_sample(W, 0, x, y)); assert(wpp >= 0.0); } assert(! isnan(wmm)); assert(! isnan(wmp)); assert(! isnan(wpm)); assert(! isnan(wpp)); /* Get the weights of the edges incident to {(x,y)}: */ double wxm = wmm + wmp; double wxp = wpm + wpp; double wym = wmm + wpm; double wyp = wmp + wpp; /* Compute the edge slopes times the edge weights: */ double wdxm = 0.0, wdxp = 0.0, wdym = 0.0, wdyp = 0.0; double wtot = wxm + wxp + wym + wyp; if (wtot == 0.0) { /* This can only happen if we asked for the full system: */ assert(full); /* Height {z[x,y]} is loose. Make it be the average of its neighbors: */ wxm = (x <= 0 ? 0.0 : 1.0); wxp = (x >= NX_G ? 0.0 : 1.0); wym = (y <= 0 ? 0.0 : 1.0); wyp = (y >= NY_G ? 0.0 : 1.0); wtot = wxm + wxp + wym + wyp; /* Leave {wdxm,wdxp,wdym,wdyp} all zero. */ } else { /* Compute {wdxm,wdxp,wdym,wdyp} from given slopes: */ wdxm = wmm*dxmm + wmp*dxmp; wdxp = wpm*dxpm + wpp*dxpp; wdym = wmm*dymm + wpm*dypm; wdyp = wmp*dymp + wpp*dypp; } assert(wtot > 0.0); assert(! isnan(wdxm)); assert(! isnan(wdxp)); assert(! isnan(wdym)); assert(! isnan(wdyp)); /* Get the index {k} into {h} of {IZ[x,y]} and its cardinal neighbors: */ long int kxm = ind_h(x-1, y); long int kxp = ind_h(x+1, y); long int kym = ind_h(x, y-1); long int kyp = ind_h(x, y+1); /* Assemble the equation for the height {Z[x,y]}: */ imgsys_equation_t *eqk = &(S->eq[k]); int nt = 0; /* Number of terms in equation. */ eqk->rhs = 0.0; eqk->ix[nt] = k; eqk->cf[nt] = 1.00; nt++; if ((wxm != 0.0) && (kxm >= 0)) { eqk->ix[nt] = kxm; eqk->cf[nt] = -wxm/wtot; eqk->rhs += +wdxm/wtot; nt++; } if ((wxp != 0.0) && (kxp >= 0)) { eqk->ix[nt] = kxp; eqk->cf[nt] = -wxp/wtot; eqk->rhs += -wdxp/wtot; nt++; } if ((wym != 0.0) && (kym >= 0)) { eqk->ix[nt] = kym; eqk->cf[nt] = -wym/wtot; eqk->rhs += +wdym/wtot; nt++; } if ((wyp != 0.0) && (kyp >= 0)) { eqk->ix[nt] = kyp; eqk->cf[nt] = -wyp/wtot; eqk->rhs += -wdyp/wtot; nt++; } eqk->nt = nt; } return S; } void pst_compute_edge_slope_and_weight ( float_image_t *G, float_image_t *GW, int axis, int x,int y, double *d, double *w ) { if (axis == 0) { /* X derivative */ pst_interpolate_four_samples(G,GW,0,x,y-2,x,y-1,x,y,x,y+1,d,w); } else if (axis == 1) { /* Y derivative */ pst_interpolate_four_samples(G,GW,1,x-2,y,x-1,y,x,y,x+1,y,d,w); } else { demand(FALSE, "invalid axis"); } return; } /* ---------------------------------------------------------------- */ /* Superseded by {float_image_mscale_mask_shrink}: */ float_image_t *pst_weight_map_shrink(float_image_t *IW, bool_t harmonic, int avgWidth) { int NC = IW->sz[0]; int NXI = IW->sz[1]; int NXJ = (NXI+1)/2; int NYI = IW->sz[2]; int NYJ = (NYI+1)/2; float_image_t *JW = float_image_new(NC, NXJ, NYJ); int xJ, yJ, c; for(yJ = 0; yJ < NYJ; yJ++) { for(xJ = 0; xJ < NXJ; xJ++) { int xI = 2*xJ, yI = 2*yJ; for (c = 0; c < NC; c++) { /* Samples in original image (with border padding): */ double w00 = get_sample(IW, c, xI, yI); double w01 = ((xI+1 < NXI) ? get_sample(IW, c, xI+1, yI) : w00); double w10 = ((yI+1 < NYI) ? get_sample(IW, c, xI, yI+1) : w00); double w11 = ((xI+1 < NXI) && (yI+1 < NYI) ? get_sample(IW, c, xI+1, yI+1) : w00); assert(w00 >= 0); assert(w01 >= 0); assert(w10 >= 0); assert(w11 >= 0); /* Average the weights: */ float w; if (harmonic) { w = 4/(1/w00 + 1/w01 + 1/w10 + 1/w11); /* Harmonic mean. */ } else { w = (w00 + w01 + w10 + w11)/4; /* Arithmetic mean. */ } assert(w >= 0); assert(! isnan(w)); /* Set sample in shrunken image: */ set_sample(JW, c, xJ, yJ, w); } } } return JW; } /* ---------------------------------------------------------------- */ /* From {pst_Kodak_Q13.c}: */ void pst_Kodak_Q13_pixel_average ( float_image_t *img, /* Image to sample. */ double xlo, /* Min X of rectangle (pixels). */ double xhi, /* Max X of rectangle (pixels). */ double ylo, /* Min Y of rectangle (pixels). */ double yhi, /* Max Y of rectangle (pixels). */ r3x3_t *P, /* Projective map from chart coords to image coords. */ double val[] /* (OUT) Pixel average. */ ); /* Extracts the average sample value from channel {c} of {img} over the axis-aligned rectangle {[xlo_xhi] × [ylo_yhi]} of the chart reference system. */ r2_t pst_Kodak_Q13_map_xy(double x, double y, r3x3_t *P); /* Applies the projective map {P} to the chart coordinates {(x,y)}, and returns the corresponding image coordinates, packaged as an {r2_t}. */ void pst_Kodak_Q13_copy_rectangle ( float_image_t *img, r3x3_t *P, double xlo, double xhi, double ylo, double yhi, int x0, int y0, int SNX, int SNY, float_image_t *DST ); void pst_Kodak_Q13_copy_rectangle ( float_image_t *img, r3x3_t *P, double xlo, double xhi, double ylo, double yhi, int x0, int y0, int SNX, int SNY, float_image_t *DST ) { if ((SNX == 0) || (SNY == 0)) { /* Nothing to do: */ return; } /* Get input image dimensions (samples): */ int NC, NX, NY; float_image_get_size(img, &NC, &NX, &NY); /* Actual pixel sizes (mm): */ double dx = (xhi - xlo)/SNX; double dy = (yhi - ylo)/SNY; /* Extract samples and store in {DST}: */ double val[NC]; int x, y; for (y = 0; y < SNY; y++) { double pix_ylo = ylo + y*dy; double pix_yhi = ylo + (y+1)*dy; for (x = 0; x < SNX; x++) { double pix_xlo = xlo + x*dx; double pix_xhi = xlo + (x+1)*dx; pst_Kodak_Q13_pixel_average(img, pix_xlo, pix_xhi, pix_ylo, pix_yhi, P, val); int c; for (c = 0; c < NC; c++) { float_image_set_sample(DST, c, x0 + x, y0 + y, val[c]); } } } } #define MAX_SUBSAMPLES (100) /* Maximum number of subsamples allowed for each output pixel. */ void pst_Kodak_Q13_pixel_average ( float_image_t *img, /* Image to sample. */ double xlo, /* Min X of rectangle (mm). */ double xhi, /* Max X of rectangle (mm). */ double ylo, /* Min Y of rectangle (mm). */ double yhi, /* Max Y of rectangle (mm). */ r3x3_t *P, /* Projective map from chart coords to image coords. */ double val[] ) { bool_t debug = FALSE; int NC = img->sz[0]; if (debug) { fprintf(stderr, " averaging"); fprintf(stderr, " [ %+6.1f _ %+6.1f ]", xlo, xhi); fprintf(stderr, " × [ %+6.1f _ %+6.1f ]\n", ylo, yhi); } /* Get corners of sampling area in {img} (homogeneous coords): */ r2_t bl = pst_Kodak_Q13_map_xy(xlo, ylo, P); r2_t br = pst_Kodak_Q13_map_xy(xhi, ylo, P); r2_t tl = pst_Kodak_Q13_map_xy(xlo, yhi, P); r2_t tr = pst_Kodak_Q13_map_xy(xhi, yhi, P); /* Get approx extent in each direction (pixels): */ double dx = fmax(r2_dist(&bl, &br), r2_dist(&tl, &tr)); double dy = fmax(r2_dist(&bl, &tl), r2_dist(&br, &tr)); /* Compute appropriate number of samples: */ double maxSamplingStep = 0.5; int nx = (int)ceil(fmax(1.0,dx)/maxSamplingStep); int ny = (int)ceil(fmax(1.0,dy)/maxSamplingStep); demand((nx <= MAX_SUBSAMPLES) && (dy <= MAX_SUBSAMPLES), "input-to-output scale is too big"); /* Compute subsampling step size (mm): */ double xStep = (xhi - xlo)/nx; double yStep = (yhi - ylo)/ny; /* Subsample and get average: */ double sumW = 0.0; double sumWV[NC], sumV2W[NC]; int c; for (c = 0; c < NC; c++) { sumWV[c] = sumV2W[c] = 0; } int ix, iy; for (iy = 0; iy < ny; iy++) { double y = ylo + (iy + 0.5)*yStep; for (ix = 0; ix < nx; ix++) { double x = xlo + (ix + 0.5)*xStep; r2_t p = pst_Kodak_Q13_map_xy(x, y, P); double W = 1.0; /* Weight of subsample. */ for (c = 0; c < NC; c++) { double V = float_image_interpolate_samples(img, c, p.c[0], p.c[1]); sumWV[c] += V*W; sumV2W[c] += V*V*W; } sumW += W; } } for (c = 0; c < NC; c++) { double avg = sumWV[c]/sumW; double var = sumV2W[c]/sumW - avg*avg; if (debug) { double dev = sqrt(fmax(var,0)); fprintf(stderr, " channel %d: avg = %8.5f dev = %8.5f\n", c, avg, dev); } val[c] = avg; } } r2_t pst_Kodak_Q13_map_xy(double x, double y, r3x3_t *P) { r3_t p = (r3_t){{ 1.0, x, y }}; r3x3_map_row(&p, P, &p); return (r2_t){{ p.c[1]/p.c[0], p.c[2]/p.c[0] }}; } /* ---------------------------------------------------------------- */ /* From {pst_sphere_basic.c}: */ /* Get the optical center {Q}: */ double Ox = O->hx; demand(! isnan(Ox), "{O.x} is NAN"); double Oy = O->hy; demand(! isnan(Oy), "{O.y} is NAN"); r2_t Q = (r2_t) {{ Ox/Om, Oy/Om }}; /* Get the displacement {Dsp} from {Q} to {K}: */ r2_t Dsp; r2_sub(S->K, &Q, &Dsp); /* Compute {D}, the distance from {Q} to {K}: */ double D = r2_norm(&Dsp); /* Compute the stretch vector {*str}: */ r2_scale((sqrt(alfa2)- 1)*radT/D, &Dsp, str); void pst_sphere_compute_ctr_rad_str_from_S_O ( pst_sphere_t *S, /* The sphere. */ hr3_point_t *O, /* Camera's viepoint. */ r2_t *ctr, /* (OUT) The center of the sphere's projection on the image. */ double *rad, /* (OUT) The minor semidiameter of that projection. */ r2_t *str /* (OUT) The stretch vector of the projection. */ ) { demand(S->R > 0, "bad sphere radius"); double Om = O->hm; demand((! isnan(Om)) && (Om >= 0), "bad {O.m}"); if (Om == 0) { /* Parallel projection (a rather common case): */ if (rad != NULL) { (*rad) = S->R; } if (ctr != NULL) { (*ctr) = S->K; } if (str != NULL) { (*str) = (r2_t) {{ 0, 0 }}; } return; } double Oz = O->hz; demand(! isnan(Oz), "{O.z} is NAN"); demand(Oz > 0, "bad {O.z}"); /* Get the center and spread of the camera: */ double G = Om/Oz; /* Compute {beta2}, the minor radius over the radius of the sphere, squared: */ double t = S->R*G, t2 = t*t; demand(t < 1, "camera is lower than top of sphere"); double beta2 = 1/(1 - t2); /* assert(beta2 >= 1); */ /* Compute the apparent transverse radius (minor semidiameter): */ double radT = S->R*sqrt(beta2); /* Worth computing, for {*rad} and {*str} */ if (radT != NULL) { (*rad) = radT; } if ((ctr != NULL) || (str != NULL)) { /* Get the optical center {Q}: */ double Ox = O->hx; demand(! isnan(Ox), "{O.x} is NAN"); double Oy = O->hy; demand(! isnan(Oy), "{O.y} is NAN"); r2_t Q = (r2_t) {{ Ox/Om, Oy/Om }}; /* Get the displacement {Dsp} from {Q} to {K}: */ r2_t Dsp; r2_sub(S->K, &Q, &Dsp); if (ctr != NULL) { /* Compute the apparent center {*ctr}: */ /* Not worth testing for {beta2 == 1}: */ r2_mix(1.0, &Q, beta2, &Dsp, ctr); } if (str != NULL) { /* Compute {D}, the distance from {Q} to {K}: */ double D = r2_norm(&Dsp); /* Compute {alfa2}, the major radius over the minor radius, squared: */ double q = D*G, q2 = q*q; double alfa2 = 1 + beta2*q2; /* Compute the stretch vector {*str}: */ /* Not worth testing for {alfa2 == 1}: */ r2_scale((sqrt(alfa2)- 1)*radT/D, &Dsp, str); } } } double *c, /* (OUT) Co-sine of elevation of {O} seen from {K}. */ double *s /* (OUT) Sine of elevation of {O} seen from {K}. */ else if (dst == +INF) { /* Oblique view from infinite distance: */ if (G != NULL) { (*G) = 0; } if (D != NULL) { (*D) = 0; } if (R != NULL) { (*R) = rad; } } /* ---------------------------------------------------------------- */ /* From {pst_sphere.c}: */ void pst_camera_args_adjust_print_viewpoint ( FILE *wr, pst_camera_t *C, double OAdj, char *fmtp, char *fmtd ) { fprintf(wr, " viewpoint"); double Oc[4] = C->O.c.c; /* Check whether it is all NAN: */ bool_t all_NAN = TRUE; int i; for (i = 0; (i < 4) && all_NAN; i++) { if (! isnan(Oc[i])) { all_NAN = FALSE; } } if (! all_NAN) { /* Must print the coords. */ /* Rescale the homogeneous coords to unit max abs: */ double amax = 0.0; for (i = 0; i < 4; i++) { double Oi = Oc[i]; if (! isnan(Oi)) { amax = fmax(amax, fabs(Oi)); } } if (amax > 0) { /* Rescale the homogeneous coords for legibility: */ hr3_point_t OR; if (! isnan(m)) { assert(isfinite(m)); double x = C->O.c.c[1]; double y = C->O.c.c[2]; double z = C->O.c.c[3]; if (m == 0) { /* Normalize point at infinity: */ double d = hypot(x, hypot(y, z)); if ((isnan(d)) || (d == 0)) { /* Give up: */ O = C->O; } else { /* Normalize {x,y,z} to unit length: */ O = (hr3_point_t) {{{ 0, x/d, y/d, z/d }}}; } } else { /* Normalize finite point to unit weight: */ O = (hr3_point_t) {{{ 1, x/m, y/m, z/m }}}; } } /* Print cleaned-up viewpoint: */ for (i = 0; i < 4; i++) { fprintf(wr, " "); fprintf(wr, (m == 0 ? fmtd : fmtp), O.c.c[i]); } } if ((! isnan(OAdj)) && (OAdj != 0)) { fprintf(wr, " adjust "); fprintf(wr, (m == 0 ? fmtd : fmtp), OAdj); } fprintf(wr, "\n"); fflush(wr); } /* ---------------------------------------------------------------- */ /* From {pst_sphere.c}: */ /* Also computes {*c = cos(e), *s = sin(e)} where {e} is the elevation angle of the viewpoint {O} as seen from the sphere's center {K}. */ void pst_sphere_compute_K_R_str_from_Q_G_dsp_rad ( r2_t *Q, /* Finite optical center on image plane. */ double G, /* Camera resolution {1/F}. */ r2_t *dsp, /* Vector from {Q} to center of ellipse. */ double rad, /* Length of minor semidiameter of ellipse. */ double *R, /* (OUT) The sphere's radius. */ r2_t *K, /* (OUT) The nominal center of the sphere. */ r2_t *str /* (OUT) Expected stretch vector of sphere. */ ) { demand(isfinite(Q->X), "invalid {Q.X}"); demand(isfinite(Q->Y), "invalid {Q.Y}"); { /* Get the direction {dir} and distance {dst} from {Q} to {ctr}: */ r2_t dir; double dst = r2_dir(&dsp, &dir); if (dst == 0) { /* Sphere on optical axis: */ (*str) = (r2_t) {{ 0, 0 }}; } else /* Compute the stretch vector {*str}: */ /* Not worth testing for {alfa2 == 1}: */ r2_scale((sqrt(alfa2)- 1)*rad, &dir, str); } } } if (debug) { if (K != NULL) { Pr(Er, " K = "); r2_print(Er, K); Pr(Er, "\n"); } if (R != NULL) { Pr(Er, " R = %18.12f\n", (*R)); } if (str != NULL) { Pr(Er, " str = "); r2_print(Er, str); Pr(Er, "\n"); } } } void pst_sphere_compute_R_dst_from_G_rad_maj (G,E->rad,maj,&(S-R),&D); void pst_sphere_compute_R_K_G_from_Q_dsp_rad_maj ( r2_t *QC, /* Finite optical center on image plane. */ r2_t *ctr, /* Center of projection. */ double rad, /* Length of minor semidiameter. */ double maj, /* Length of major semidiameter. */ double *R, /* (OUT) The sphere's radius. */ r2_t *K, /* (OUT) The nominal center of the sphere. */ double *G /* (OUT) Camera resolution inferred from {E,Q}. */ ) { bool_t debug = TRUE; demand(E->rad > 0); demand(isfinite(QC->X), "invalid {QC.X}"); demand(isfinite(QC->Y), "invalid {QC.Y}"); /* Compute the displacement {dsp} from {QC} to {ctr}: */ r2_t dsp; r2_sub(ctr, &Q, &dsp); double G = 1/O->Z; /* Compute the the ratio {u=rad/F} where {F} is the focal length: */ double u = rad*G, u2 = u*u; /* Compute {beta2}, the minor radius over the sphere's radius, squared: */ double beta2 = 1 + u2; if (debug) { Pr(Er, " beta2 = %16.12f beta = %16.12f\n", beta2, sqrt(beta2)); } /* assert(beta2 >= 1); */ /* Compute the sphere's radius {*R}: */ if (R != NULL) { (*R) = rad/sqrt(beta2); } if ((K != NULL) || (str != NULL)) { /* Get the displacement {dsp} from {Q} to {ctr}: */ if (K != NULL) { /* Compute the nominal center {*K}: */ /* Not worth testing for {beta2 == 1}: */ r2_mix(1.0, &Q, 1/beta2, &dsp, K); } if (str != NULL) { /* Get the direction {dir} and distance {dst} from {Q} to {ctr}: */ r2_t dir; double dst = r2_dir(&dsp, &dir); if (dst == 0) { /* Sphere on optical axis: */ (*str) = (r2_t) {{ 0, 0 }}; } else { /* Compute {alfa2}, the major radius over the minor radius, squared: */ double v = dst*G, v2 = v*v; double alfa2 = 1 + beta2*v2; if (debug) { Pr(Er, " alfa2 = %16.12f alfa = %16.12f\n", alfa2, sqrt(alfa2)); } /* Compute the stretch vector {*str}: */ /* Not worth testing for {alfa2 == 1}: */ r2_scale((sqrt(alfa2)- 1)*rad, &dir, str); } } } if (debug) { if (K != NULL) { Pr(Er, " K = "); r2_print(Er, K); Pr(Er, "\n"); } if (R != NULL) { Pr(Er, " R = %18.12f\n", (*R)); } if (str != NULL) { Pr(Er, " str = "); r2_print(Er, str); Pr(Er, "\n"); } } } /* ---------------------------------------------------------------- */ /* From {pst_geom.c}: */ ellipse_crs_t *pst_geom_sphere_new(void); /* Alocates a new {elipse_t} record, not initialized. */ ellipse_crs_t *pst_geom_sphere_new(void) { ellipse_crs_t *E = (ellipse_crs_t *)notnull(malloc(sizeof(ellipse_crs_t)), "no mem"); E->ctr = (r2_t){{ 0, 0 }}; E->rad = 0; E->str = (r2_t){{ 0, 0 }}; return E; } r3_t pst_geom_sphere_compute_normal(r2_t *uv); /* Computes the outwards normal vector to the sphere's surface at the point that projects onto the point {uv}. Assumes a {U,V,W} coordinate system where the sphere has radius 1 and center at the origin, and {W} points towards the camera. !!! Use true perspective instead of affine approximation. !!! */ void pst_geom_sphere_view_matrices ( ellipse_crs_t *E, /* Geometry of sphere's projection. */ r3x3_t *XY1_to_UV1, r3x3_t *UVW_to_XYZ ); /* !!! Use true perspective instead of affine approximation. !!! Computes the projection matrices {xym_to_uvm} and {uvw_to_xyz} needed by {pst_normal_map_from_proc}. Besides the Image Coordinate System (ICS), these procedures use another orthogonal three-dimensional coordinate system called the Object Coordinate System (OCS). The OCS coordinates {(U,V,W)} are measured in multiples of the sphere's radius. The origin {(0,0,0)} is the the sphere's center. The {U} and {V} axes are perpendicular to the viewing direction, so that their tips are silhouette points --- that is, their projections lie on the outline {E} of the sphere's projection. The {U} axis is parallel to the image plane, so that the points {(1,0,0)} and {(0,1,0)} of the OCS project onto the tips of a minor and a major semidiameter of {E}, respectively; and {(0,0,1)} is the point on the sphere's surface that is closest to the camera, that, like {(0,0,0)}, projects onto the center of {E}. It follows that, in OCS coordinates, the sphere in question has unit radius and is centered at the origin. Thus, the {W} coordinate of the visible point of the sphere with given {U,V} coordinates is {W = sqrt(1-U^1-V^2)} that is, {W = pst_geom_sphere_compute_normal((U,V)}. Moreover, sphere's unit normal vector at a point with OCS coordinates {(U,V,W)} is {(U,V,W)}. The matrix {xym_to_uvm} computed by this procedure is an affine map that takes the homogeneous ICS coordinates {[X,Y,1]} of a point of the image to the homogneous OCS coordinates {[U,V,1]} of a point on the {W=0} plane that projects onto the image point {(X,Y)}. The second matrix, {uvw_to_xyz}, is a linear map that converts Cartesian OCS coordinates {(U,V,W)} of a OCS unit vector to the ICS coordinates {(X,Y,Z)} of the ICS unit vector with same direction. (Note that this is a /linear/ map for /vectors/, not an /affine/ map for /points/.) Thus, to obtain the true normal at the point with image coordinates {(X,Y)}, one should use {xym_to_uvm} to compute the corresponding {(U,V)} coordinates, then {pst_geom_sphere_compute_normal} to compute the {W} coordinate, then {uvw_to_xyz} to get the normal's coordinates in the ICS. */ void pst_geom_sphere_view_matrices ( ellipse_crs_t *E, /* Geometry of sphere's projection. */ r3x3_t *xym_to_uvm, r3x3_t *uvw_to_xyz ) { bool_t debug = TRUE; /* Initialize {xym_to_uvm} with translation matrix that subtracts {E->ctr}: */ r3x3_ident(xym_to_uvm); xym_to_uvm->c[2][0] = -E->ctr.c[0]; xym_to_uvm->c[2][1] = -E->ctr.c[1]; /* Initialize {uvw_to_xyz} with the identity: */ r3x3_ident(uvw_to_xyz); /* Compute min and max radii: */ double mstr = r2_norm(&(E->str)); double rmin = E->rad; double rmax = E->rad + mstr; if (rmin != rmax) { /* Elliptic shape, must change basis: */ /* Get directions {du,dv} of shortest and longest axes: */ r2_t dmin, dmax; (void)r2_dir(&(E->str), &dmax); dmin = (r2_t){{ +dmax.c[1], -dmax.c[0] }}; /* Append to {xym_to_uvm} a rotation that maps {x,y} to {u,v} */ r3x3_t rot; r3x3_ident(&rot); rot.c[0][0] = dmin.c[0]; rot.c[1][0] = dmin.c[1]; rot.c[0][1] = dmax.c[0]; rot.c[1][1] = dmax.c[1]; r3x3_mul(xym_to_uvm, &rot, xym_to_uvm); /* Set {uvw_to_urz} to rotate around {u} axis, changing {u,v,w} into {u,r,z}: */ r3x3_t uvw_to_urz; r3x3_ident(&uvw_to_urz); double ct = rmin/rmax; /* Cosine of rotation angle. */ double st = sqrt(1.0 - ct*ct); /* Sine of rotation angle. */ uvw_to_urz.c[1][1] = +ct; uvw_to_urz.c[1][2] = +st; uvw_to_urz.c[2][1] = -st; uvw_to_urz.c[2][2] = +ct; /* Append to {uvw_to_urz} the inverse rotation from {u,r,z} to {x,y,z}: */ r3x3_t tor; r3x3_ident(&tor); tor.c[0][0] = dmin.c[0]; tor.c[0][1] = dmin.c[1]; tor.c[1][0] = dmax.c[0]; tor.c[1][1] = dmax.c[1]; r3x3_mul(&uvw_to_urz, &tor, uvw_to_xyz); } /* Append to {xym_to_uvm} an unscale by {rmin,rmax}: */ r3x3_t unrad; r3x3_ident(&unrad); unrad.c[0][0] = 1/rmin; unrad.c[1][1] = 1/rmax; r3x3_mul(xym_to_uvm, &unrad, xym_to_uvm); if (debug) { fprintf(stderr, "xym_to_uvm:\n"); r3x3_print(stderr, xym_to_uvm); fprintf(stderr, "\n"); fprintf(stderr, "uvw_to_xyz:\n"); r3x3_print(stderr, uvw_to_xyz); fprintf(stderr, "\n"); } } r3_t pst_geom_sphere_compute_normal(r2_t *uv) { double u = uv->c[0]; double v = uv->c[1]; /* Check if inside unit {u,v} disk: */ double uv2 = u*u + v*v; if (uv2 >= 0.999999) { /* Outside - return null normal: */ return (r3_t) {{ 0.0, 0.0, 0.0 }}; } else { /* Inside - compute {w} coordinate: */ double w = sqrt(1.000001 - uv2); /* The normal is the points's position: */ return (r3_t) {{ u, v, w }}; } } /* ---------------------------------------------------------------- */ /* From {pst_geom.c}: */ double pst_camera_distance_from_resolution(double res); /* Computes the distance {z} from the camera to the `true' image plane (where pixels are 1 unit wide), given the angular resolution {res} at the optical axis. */ double pst_camera_resolution_from_distance(double z); /* Computes the camera's angular resolution {res} at the optical axis given its distance from the camera to the `true' image plane (where pixels are 1 unit wide). It is the inverse of {pst_camera_distance_from_resolution}. */ double pst_camera_distance_from_resolution(double res) { return 2/tan(0.5*fabs(res)); } double pst_camera_resolution_from_distance(double z); { return 2*atan(0.5/fabs(z)); } /* ---------------------------------------------------------------- */ /* From {pst_geom.c}: */ void pst_geom_sphere_compute_stretch_from_apparent_params ( r2_t *opa, /* Optical axis of camera. */ double res, /* Angular resolution of camera. */ r2_t *ctr, /* The center of the sphere's projection on the image. */ double rad, /* The minor semidiameter of the projection. */ r2_t *strP /* (OUT) The major semidiameter of the projection. */ ) { r2_t tct; /* The `true' center. */ double trd; /* The `true' radius. */ pst_geom_sphere_compute_true_params_from_apparent_params ( opa, res, ctr, rad, &tct, &trd ); pst_geom_sphere_compute_stretch_from_true_params ( opa, res, &tct, trd, strP ); } /* ---------------------------------------------------------------- */ /* From {pst_camera.c}: */ /* Try to determine the natural width of {fmtp,fmtr}-printed values: */ char *tmp = NULL; asprintf(&tmp, fmtp, 0); int wdp = strlen(tmp); free(tmp); tmp = NULL; asprintf(&tmp, fmtr, 0); int wdr = strlen(tmp); free(tmp); /* ---------------------------------------------------------------- */ /* From {pst_geom.h}: */ /* CAMERA PARAMETERS */ #define pst_geom_aperture_HELP \ "-aperture {APERTURE}" #define pst_geom_aperture_HELP_INFO \ "This option specifies the aperture of the camera," \ " defined as {tan(THETA/2)} where {THETA} is the" \ " horizontal angular width of the camera's" \ " field of view. Thus, if an object of width {W} occupies" \ " the whole width of the image" \ " when placed at distance {D} from the camera, the" \ " {APERTURE} is {W/D/2}." #define pst_geom_optical_axis_HELP \ "-opticalAxis {OPX} {OPY}" #define pst_geom_optical_axis_HELP_INFO \ "This option specifies the coordinates of the camera's" \ " optical axis in the image, as being {(OPX,OPY)}." /* PRINTOUT */ /* ---------------------------------------------------------------- */ /* From {pst_geom.c}: */ r2_t pst_geom_sphere_est_stretch(r2_t *ctr, double rad, int NX, int NY, double ap); /* Computes the stretch vector for a sphere whose projection has center {ctr} and minor semidiameter {rad} (both in pixels, in the image coordinate system). Assumes that the image has {NX} columns by {NY} rows, with pixel in column {x} and row {y} being the square {[x_x+1]×[y_y+1]}. Assumes also that the optical axis goes through the center of the image, and that {ap} is {tan(theta/2)} where {theta} is the total angular width of the the camera's field of view in the horizontal direction. */ r2_t pst_geom_sphere_est_stretch(r2_t *ctr, double rad, int NX, int NY, double ap) { r2_t str = (r2_t) {{ 0, 0 }}; if (ap == 0) { /* Camera at infinity: */ return str; } /* In what follows, the /image plane/ is the plane perpendicular */ /* to the camera axis that passes through the sphere's center. */ r2_t h = (r2_t){{ 0.5*NX, 0.5*NY }}; /* Half-width and center of image. */ r2_t c; r2_sub(ctr, &h, &c); /* Sphere ctr disp from image center. */ r2_t cv; double cr = r2_dir(&c, &cv); /* Direction from img ctr to sphere ctr. */ if (cr < 1.0e-100) { /* Sphere on optical axis: */ return str; } double F = h.c[0]/ap; /* Distance of camera from image plane. */ double t = cr/F; /* Tangent of view tilt angle. */ double RAD = rad*hypot(1, t); /* Major semidiameter or ellipse. */ r2_scale(RAD - rad, &cv, &str); return str; } /* ---------------------------------------------------------------- */ /* From {pst_fit_ellipse.c}: */ void pst_fit_ellipse_adjust_bounds(double *mid, double *adj, double low); /* Adjusts the arguments {*mid} and {*adj} so as to clip the range {[mid-adj _ mid+adj]} against the range {[low _ +INF]}. In particular, if {mid+adj} is less than {low}, sets {mid} to {low} and {adj} to 0. The value of {adj} must be non-negative. */ void pst_fit_ellipse_adjust_bounds(double *mid, double *adj, double low) { double md = (*mid), ad = (*adj); if (md + ad <= low) { md = low; ad = 0.0; } else if (md - ad < low) { md = (low + md + ad)/2; ad = md - low; } (*mid) = md; (*adj) = ad; } /* Adjust {radius}, {radAdj} so that trial radii are not too small: */ double radMin = 2.0; /* Minimum sphere radius in pixels. */ pst_fit_ellipse_adjust_bounds(&(locE.rad), &radAdj, radMin); /* Adjust {strLen}, {strAdj} so as to avoid negative stretch: */ double strLen = r2_norm(&(locE.str)); pst_fit_ellipse_adjust_bounds(&strLen, &strAdj, 0.0); /* ---------------------------------------------------------------- */ /* From {pst_geom.c}: */ void pst_geom_sphere_int_bbox ( ellipse_t *ep, double mrg, /* Extra margin. */ int *xLoP, /* (OUT) Min X of clip area. */ int *xHiP, /* (OUT) Max X of clip area. */ int *yLoP, /* (OUT) Min Y of clip area. */ int *yHiP /* (OUT) Max Y of clip area. */ ) { double xLoF, xHiF, yLoF, yHiF; float_image_paint_get_ellipse_bbox ( ep->ctr.c[0], ep->ctr.c[1], ep->rad, ep->str.c[0], ep->str.c[1], &xLoF, &xHiF, &yLoF, &yHiF ); (*xLoP) = imax(0, (int)floor(xLoF - mrg)); (*xHiP) = imin(nx-1, (int)ceil (xHiF + mrg)); (*yLoP) = imax(0, (int)floor(yLoF - mrg)); (*yHiP) = imin(ny-1, (int)ceil (yHiF + mrg)); } /* ---------------------------------------------------------------- */ /* From {pst_fit_ellipse.c}: */ float_image_t *pst_fit_ellipse_rasterize ( int xLo, int xHi, int yLo, int yHi, ellipse_t *ep ); /* Creates a rasterized image of the ellipse described by {ep}. The image is 1.0 inside the ellipse, 0.0 outside. Assumes that the ellipse's width and length are large compared to the pixel size. The image starts at column {xLo} and row {yLo} of the full image, and has {xHi-xLo) columns and {yHi-yLo} rows. */ float_image_t *pst_fit_ellipse_rasterize ( ellipse_t *ep, int xLo, int xHi, int yLo, int yHi ) { /* Find a box around the ellipse, with ~2 pixels extra on each side: */ int xLo, xHi, yLo, yHi; float_image_paint_get_ellipse_bbox ( xctr, yctr, erad, xstr, ystr, &xLo, &xHi, &yLo, &yHi ); xLo = imax(0, xLo - 2); xHi = imin(nx-1, xHi + 2); yLo = imax(0, yLo - 2); yHi = imin(ny-1, yHi + 2); int NX = xHi - xLo; int NY = yHi - yLo; /* Compute the ellipse's center {xctr,yctr} rel to bottom left corner: */ double xctr = ep->ctr.c[0] - xLo; double yctr = ep->ctr.c[1] - yLo; double rad = ep->rad; double xstr = ep->str.c[0]; double ystr = ep->str.c[1]; float_image_t *EMG = float_image_new(1, NX, NY); float_image_paint_ellipse ( EMG, 0, 0, NX-1, 0, NY-1, xctr, yctr, rad, xstr, ystr, 0.0, 1.0, 2 ); return EMG; } /* ---------------------------------------------------------------- */ /* From {pst_fit_ellipse.c}: */ double pst_fit_ellipse_eval(float_image_t *IGR, pst_geom_sphere_t *geo) { /* Get the dimensions of the gradient image: */ demand(IGR->sz[0] == 1, "image must be monochromatic"); int NX = IGR->sz[1]; int NY = IGR->sz[2]; /* Rasterize an image of the solid ellipse: */ int xLo, xHi, yLo, yHi; float_image_paint_get_ellipse_bbox ( float_image_t *EMG = pst_fit_ellipse_rasterize ( geo, &xLo, &xHi, &yLo, &yHi ); /* Compute its relative gradient: */ float_image_t *EGR = float_image_gradient_sqr_relative(EMG, 0.01); /* Compute the weight mask: */ int hw = 4; assert(2*hw+1 == 9); double D = 70; double wd[9] = { 1/D, 8/D, 28/D, 56/D, 70/D, 56/D, 28/D, 1/D }; float_image_t *EMK = float_image_mmorph_dilate(EGR, hw, wd); /* Scan the box and accumulate the differences squared: */ double sum_w = 0; double sum_w_ds2 = 0; int xI, yI; for (xI = xLo; xI < xHi; xI++) { int xE = xI - xLo; for (yI = yLo; yI < yHi; yI++) { int yE = yI - yLo; /* Get the value {sI} of the pixel in the {IGR} image: */ double sI = float_image_get_sample(IGR, 0, xI, yI); /* Get the value {sE} of the pixel in the {EGR} image: */ double sE = float_image_get_sample(EGR, 0, xE, yE); /* Get the weight {w} of the pixel in the ideal grad map: */ double w = float_image_get_sample(EMK, 0, xE, yE); if (w > 0.0) { double ds = sI - sE; sum_w += w; sum_w_ds2 += w*ds*ds; } } } float_image_free(EMK); float_image_free(EGR); float_image_free(EMG); assert(sum_w > 0); return sum_w_ds2 / sum_w; } /* ---------------------------------------------------------------- */ /* From {pst_gray_scale_fit.c}: */ /* Get corners of sampling area in {img} (homogeneous coords): */ r2_t bl = pst_gray_scale_fit_map_xy(cX-rx, cY-ry, P); r2_t br = pst_gray_scale_fit_map_xy(cX+rx, cY-ry, P); r2_t tl = pst_gray_scale_fit_map_xy(cX-rx, cY+ry, P); r2_t tr = pst_gray_scale_fit_map_xy(cX+rx, cY+ry, P); /* The following code assumes that the convex quadrilateral {Q} with corners {bl,br,tl,tr} in {img} is not too far from a rectangle, so that a uniform rectangular grid in chart space provides adequate sampling in image space. If {Q} is a highly sheared paralleogram, the procedure will waste more samples than necessary. If {Q} has very unequal opposite sides, the density of sampling points will be non-uniform. In either case, the only consequences will be wasted effort and/or a slightly sub-optimal sensitivity to image noise. */ /* Get approx extent in each direction (pixels): */ double dx = fmax(r2_dist(&bl, &br), r2_dist(&tl, &tr)); double dy = fmax(r2_dist(&bl, &tl), r2_dist(&br, &tr)); /* Compute appropriate number of samples on each side of center: */ double maxStep = 0.5; int nx = (int)ceil(fmax(1.0,dx/2)/maxStep); int ny = (int)ceil(fmax(1.0,dy/2)/maxStep); /* Compute subsampling step size (mm): */ double xStep = C.rx/nx; double yStep = C.ry/ny; /* ---------------------------------------------------------------- */ /* From {pst_Kodak_Q13.h}: */ #define pst_Kodak_Q13_scale_width \ (2*pst_Kodak_Q13_end_patch_width + (pst_Kodak_Q13_num_steps-2)*pst_Kodak_Q13_mid_patch_width) /* Total X extent of gray scale, including the first and last patches (mm). */ #define pst_Kodak_Q13_scale_start_x \ ((pst_Kodak_Q13_total_width - pst_Kodak_Q13_scale_width)/2) /* X position of left edge of first gray patch (mm). */ /* Help/info for radius-only lamp specification: */ /* ---------------------------------------------------------------- */ /* THE INTEGRATION SYSTEM The procedures in this section store into the {imgsys_t} system {S} the linear equations that yield the height map {IZ} given a slope map {IDZ}. The value of{IDZ[0,x,y]} must be the X derivative {dZ/dX} of the heigh map, averaged over the pixel {[x,y]} The value of {IDZ[1,x,y]} must be the Y derivative {dZ/dY} averaged over the same pixel. Some of these procedures also take a weight map {IW}, with the same number of rows and columns as the slope map {IDZ}. The value of {IW[0,x,y]} is interpreted as the reliability of the slope data {IDZ[c,x,y]}. The absolute value of the weights is not relevant, only their ratios. In particular, if {IW[0,x,y]} is zero, the slope in pixel {[x,y]} is considered indeterminate. The system has the form {M h = b} where {M} is a known matrix, {h} is the unknown solution vector, and {b} is a known vector. The unknowns are the heights (Z coordinates) at the CORNERS of all pixels in {IDZ}. Therefore, if the slope map {IDZ} has {NX} columns and {NY} rows, the system has {NH = (NX+1)*(NY+1)} equations and unknowns. More precisely, {h[x + (NX+1)*y} is the height at the integer point {(x,y)}, not at the pixel center {(x+0.5,y+0.5)}, for all {x} in {0..NX} and all {y} in {0..NY}. */ void pst_slope_map_build_integration_system(float_image_t *IDZ, imgsys_t* S); /* Stores into {S} the slope integration system for the given slope map {IDZ} . The equations state that the Laplacian computed from the (unknown) height map {IZ} is equal to that computed from the (known) slope maps {IDZ}. Different equations, based on equality of slopes, are used along the edges. The system has a one-dimensional family of solutions. The homogeneous system {M h = 0} has a nonzero solution, namely {h[k]==1} for all {k}. */ void pst_slope_map_build_integration_system_1(float_image_t *IDZ, imgsys_t* S) { /* Each equation and each unknown of {S} corresponds to one element in the height map {IZ}, that is, to the {Z} coordinate at a vertex (pixel corner) of the slope map {IDZ}. By definition, the vertex {IZ[x,y]} has coordinates {(x,y)}. Therefore the pixel with samples {IDZ[c,x,y]} (for {c=0,1}) is the unit square whose lower left corner is the vertex with indices {[x,y]}, that is, whose diagonal goes from {(x,y)} to {(x+1,y+1)}. In general, the equation of {S} associated to the vertex {[x,y]} has the form {Lz[x,y] = Ld[x,y]}, where {Lz,Ld} are two estimates of the Laplacian {d2z/dx2 + d2z/dy2} of the height function {z(x,y)}, evaluated at the vertex {(x,y)}. The {Lz} estimate is obtained from the height map {IZ}, by finite differences between the height {IZ[x,y]} at that vertex and the heights of its four nearest neighbors; and therefore is a linear function of those five unknowns. The {Ld} estimate is obtained from the given slope map {IDZ}, evaluated at the centers of the four pixels incident to the vertex {[x,y]}; that is, from the entries {IDZ[c,x-dx,y-dy]}, where {dy}, {dx}, and {c} are 0 or 1. Observe that {Ld} is independent from the unknowns. Actually, this description holds only for the vertices {[x,y]} that are surrounded by four pixels in the slope mape {IDZ} For the vertices with indices {[x,0]}, along the lower edge of the images (but excluding the corners {[0,0]} and {[NX,0]}), we use the equation {Dz[x,0] = Dd[x,0]}, where {Dz} and {Dd} are two estimates of the {Y} derivative of the height function at the point {(x,0.5)}. The {Dz} estimate is obtained by finite differences between the (unknown) heights of the two nearest vertices, that is, {Dz[x,0] = IZ[x,1] - IZ[x,0]}. The {Dd} estimate is obtained from the given slopes {IDZ[c,x-dx,y]} where {dx} is 0 or 1. The equations along the other three edges are similar. At the vertex {[0,0]} (lower left corner of the slope map {IDZ}), the equation has the same form {Dz[x,y] = Dd[x,y]}, except that the height estimates {Dz} and {Dd} apply to the point {(0.5,0.5)}, and the derivative is taken with respect to the vector {(1,1)}. That is, {Dz[0,0] = IZ[1,1] - IZ[0,0]}, and {Dd[0,0] = IDZ[0,0,0] + IDZ[1,0,0]}. The equations for the other three image corners (vertices {[NX,0]}, {[0,NY]} and {[NX,NY]}) are similar. This system will be solved by the Gauss (or Gauss-Jordan) iterative method. It is an underconstrained system, since the homogeneous version {M Z = 0} has a non-zero solution --- namely, all heights equal to 1. (Note that evert row of {M} adds to zero, since it is an estimator for the derivative or Laplacian of the height field.) */ /* Get/check the sizes of the slope maps: */ int NC = IDZ->sz[0]; assert(NC == 2); int NX = IDZ->sz[1]; int NY = IDZ->sz[2]; /* Check the size of the system: */ int NH = (NX+1)*(NY+1); /* Number of unknowns (heights). */ assert(S->N == NH); auto long int ind_h(int x, int y); /* This function computes the index into the solution vector {h} corresponding to vertex {[x,y]} of the images. */ long int ind_h(int x, int y) { return x + y*(NX+1); } /* Make sure that we can build the system: */ assert(MAXCOEFS >= 5); /* Enumerate the vertices of {IDX,IDY}, i.e. elems of the height map: */ int x, y; for(y = 0; y <= NY; y++) { for(x = 0; x <= NX; x++) { /* Get the index {k} into {h} of {IZ[x,y]}: */ long int k = ind_h(x, y); /* Get equation {k} and clear it out: */ imgsys_equation_t *eqi = &(S->eq[k]); int kk; for (kk = 0; kk < MAXCOEFS; kk++) { eqi->coeff[kk] = 0.0; eqi->ix[kk] = -1; } eqi->rhs = 0.0; /* Build the equation: */ if ((y == 0) && (x == 0)) { /* Corner vertex {(0,0)}: */ eqi->ix[0] = k; eqi->coeff[0] = - 1.00; eqi->ix[1] = ind_h(1,1); eqi->coeff[1] = + 1.00; float dx00 = get_sample(IDZ, 0, 0, 0); float dy00 = get_sample(IDZ, 1, 0, 0); eqi->rhs = + dx00 + dy00; } else if ((y == 0) && (x == NX)) { /* Corner vertex {(NX,0)}: */ eqi->ix[0] = k; eqi->coeff[0] = - 1.00; eqi->ix[1] = ind_h(NX-1,1); eqi->coeff[1] = + 1.00; float dx10 = get_sample(IDZ, 0, NX-1, 0); float dy10 = get_sample(IDZ, 1, NX-1, 0); eqi->rhs = - dx10 + dy10; } else if ((y == NY) && (x == 0)) { /* Corner vertex {(0,NY)}: */ eqi->ix[0] = k; eqi->coeff[0] = - 1.00; eqi->ix[1] = ind_h(1,NY-1); eqi->coeff[1] = + 1.00; float dx01 = get_sample(IDZ, 0, 0, NY-1); float dy01 = get_sample(IDZ, 1, 0, NY-1); eqi->rhs = + dx01 - dy01; } else if ((y == NY) && (x == NX)) { /* Corner vertex {(NX,NY)}: */ eqi->ix[0] = k; eqi->coeff[0] = - 1.00; eqi->ix[1] = ind_h(NX-1,NY-1); eqi->coeff[1] = + 1.00; float dx11 = get_sample(IDZ, 0, NX-1, NY-1); float dy11 = get_sample(IDZ, 1, NX-1, NY-1); eqi->rhs = - dx11 - dy11; } else if (y == 0) { /* Border vertex {(x,0)}: */ eqi->ix[0] = k; eqi->coeff[0] = - 1.00; eqi->ix[1] = ind_h(x,1); eqi->coeff[1] = + 1.00; float dy00 = get_sample(IDZ, 1, x, 0); float dy10 = get_sample(IDZ, 1, x-1, 0); eqi->rhs = + (dy00 + dy10)/2; } else if (y == NY) { /* Border vertex {(x,NY)}: */ eqi->ix[0] = k; eqi->coeff[0] = - 1.00; eqi->ix[1] = ind_h(x,NY-1); eqi->coeff[1] = + 1.00; float dy01 = get_sample(IDZ, 1, x, NY-1); float dy11 = get_sample(IDZ, 1, x-1, NY-1); eqi->rhs = - (dy01 + dy11)/2; } else if (x == 0) { /* Border vertex {(0,y)}: */ eqi->ix[0] = k; eqi->coeff[0] = - 1.00; eqi->ix[1] = ind_h(1,y); eqi->coeff[1] = + 1.00; float dx00 = get_sample(IDZ, 0, 0, y); float dx01 = get_sample(IDZ, 0, 0, y-1); eqi->rhs = + (dx00 + dx01)/2; } else if (x == NX) { /* Border vertex {(NX,y)}: */ eqi->ix[0] = k; eqi->coeff[0] = - 1.00; eqi->ix[1] = ind_h(NX-1,y); eqi->coeff[1] = + 1.00; float dx10 = get_sample(IDZ, 0, NX-1, y); float dx11 = get_sample(IDZ, 0, NX-1, y-1); eqi->rhs = - (dx10 + dx11)/2; } else { /* Interior vertex {[x,y]}: */ eqi->ix[0] = k; eqi->coeff[0] = - 4.00; eqi->ix[1] = ind_h(x,y-1); eqi->coeff[1] = + 1.00; eqi->ix[2] = ind_h(x,y+1); eqi->coeff[2] = + 1.00; eqi->ix[3] = ind_h(x-1,y); eqi->coeff[3] = + 1.00; eqi->ix[4] = ind_h(x+1,y); eqi->coeff[4] = + 1.00; float dx00 = get_sample(IDZ, 0, x, y); float dx01 = get_sample(IDZ, 0, x, y-1); float dx10 = get_sample(IDZ, 0, x-1, y); float dx11 = get_sample(IDZ, 0, x-1, y-1); float dy00 = get_sample(IDZ, 1, x, y); float dy01 = get_sample(IDZ, 1, x, y-1); float dy10 = get_sample(IDZ, 1, x-1, y); float dy11 = get_sample(IDZ, 1, x-1, y-1); eqi->rhs = (dx00 - dx10 + dx01 - dx11 + dy00 - dy01 + dy10 - dy11)/2; } } } } void pst_slope_map_build_integration_system_2(float_image_t *IDZ, float_image_t *W, imgsys_t* S); /* Stores into {S} the slope integration system for the given slope map {IDZ} . The equations state that the height difference along each pixel diagonal must be equal to the slope along that diagonal. These constraints are weighted by the slope reliability weights. The system splits into two connected Each connected component of the system has a -dimensional family of solutions. The homogeneous system {M h = 0} has a nonzero solution, namely {h[k]==1} for all {k}. */ void pst_slope_map_build_integration_system_2(float_image_t *IDZ, float_image_t *IW, imgsys_t* S) { /* Get/check the sizes of the slope maps: */ int NC = IDZ->sz[0]; assert(NC == 2); int NX = IDZ->sz[1]; int NY = IDZ->sz[2]; /* Check the size of the system: */ int NH = (NX+1)*(NY+1); /* Number of unknowns (heights). */ assert(S->N == NH); /* The unknowns are the list {h} of all the {NH} heights, namely {h[x + (NX+1)*y]$ is the height {z[x,y]} at vertex {(x,y)}. Conceptually, we define a quadratic mismatch term {q[x,y](h)} for each pixel {x,y} in the slope map {IDZ}. The term is the square of the difference between the (known) mean gradient within the pixel given by {IDZ}, and the mean gradient computed by numerical difference of the (unknown) heights at the four corners of the pixel: { q[x,y](h) = (dx - (- h[k00] - h[k01] + h[k10] + h[k11])/2)^2 + (dy - (- h[k00] + h[k01] - h[k10] + h[k11])/2)^2 } where {dx = IDZ[0,x,y]} and {dy = IDZ[1,x,y]} are the given mean slopes, and {h[k00]}, {h[k01]}, {h[k10]}, and {h[k11]} are the {h} variables corresponding to heights {z[x,y]}, {z[x,y+1]}, {z[x+1,y]}, and {z[x+1,y+1]}, respectively. We want to minimize the sum {Q(h)} of all terms {q[x,y](h)}, each weighted bt {W[x,y]}. This is equivalent to finding the {h} vector that anihilates the gradient {G} of {Q}. The gradient of {Q} is an affine function of {h}, namely {G(h) = A h - b}. Each term {W[x,y]*q[x,y](h)} of {Q} contributes to four components of the gradient, namely those with coordinates {k00,k01,k10,k11}. These contributions are: row {k00}: { + h[k00] - h[k11] - (- dx - dy) } row {k01}: { + h[k01] - h[k10] - (- dx + dy) } row {k10}: { - h[k01] + h[k10] - (+ dx - dy) } row {k11}: { - h[k00] + h[k11] - (+ dx + dy) } all multiplied by {W[x,y]}. Thus, to build the system, we add those contributions to the appropriate elements of {A} and {b}. Note that if {dx == dy == 0}, all those terms are satisfied by the two solutions (1) all heights equal to 1, and (2) all heights either +1 or -1, in checkerboard fashion. Therefore, if the system is to be solved iteratively, one must take care to exclude those two solutions. */ auto long int ind_h(int x, int y); /* This function computes the index into the solution vector {h} corresponding to vertex {[x,y]} of the images. */ long int ind_h(int x, int y) { return x + y*(NX+1); } /* Make sure that we can build the system: */ assert(MAXCOEFS >= 5); /* Enumerate the pixels of the height map: */ int x, y; for(y = 0; y <= NY; y++) { for(x = 0; x <= NX; x++) { /* Get the indices into {h} of {IZ[x,y]} and its four diagonal neighbors: */ long int k = ind_h(x, y); long int kmm = ind_h(x-1, y-1); long int kmp = ind_h(x-1, y+1); long int kpm = ind_h(x+1, y-1); long int kpp = ind_h(x+1, y+1); /* Get equation {k} and clear it out: */ imgsys_equation_t *eqi = &(S->eq[k]); int kk; for (kk = 0; kk < MAXCOEFS; kk++) { eqi->coeff[kk] = 0.0; eqi->ix[kk] = -1; } eqi->rhs = 0.0; /* Get the diagonal slopes and weights in the four adjoining pixels: */ double dmm, dmp, dpm, dpp = 0.0; double wmm, wmp, wpm, wpp = 0.0; if ((x > 0) && (y > 0)) { dmm = + get_sample(IDZ, 0, x-1, y-1) + get_sample(IDZ, 1, x-1, y-1); wmm = (IW == NULL ? 1.0 : get_sample(IW, 0, x-1, y-1)); } if ((x < NX-1) && (y > 0)) { dpm = - get_sample(IDZ, 0, x+1, y-1) + get_sample(IDZ, 1, x+1, y-1); wpm = (IW == NULL ? 1.0 : get_sample(IW, 0, x+1, y-1)); } if ((x > 0) && (y < NY-1)) { dmp = + get_sample(IDZ, 0, x-1, y+1) - get_sample(IDZ, 1, x-1, y+1); wmp = (IW == NULL ? 1.0 : get_sample(IW, 0, x-1, y+1)); } if ((x < NX-1) && (y < NY-1)) { dpp = - get_sample(IDZ, 0, x+1, y+1) - get_sample(IDZ, 1, x+1, y+1); wpp = (IW == NULL ? 1.0 : get_sample(IW, 0, x+1, y+1)); } /* Build the equation: */ eqi->ix[0] = k; eqi->coeff[0] = + wmm + wmp + wpm + wpp; eqi->ix[1] = kmm; eqi->coeff[1] = - wmm; eqi->ix[2] = kmp; eqi->coeff[2] = - wmp; eqi->ix[3] = kpm; eqi->coeff[3] = - wpm; eqi->ix[4] = kpp; eqi->coeff[4] = - wpp; eqi->rhs = wmm*dmm + wmp*dmp + wpm*dpm + wpp*dpp; } } } /* ---------------------------------------------------------------- */ /* Compute total strip image dimensions and starting row of each strip: */ int Y0_r0 = 0; int Y0_q0 = Y0_r0 + SNY_r0; int Y0_x0 = Y0_q0 + SNY_q0; int Y0_s0 = Y0_x0 + SNY_x0; int Y0_s1 = Y0_s0 + SNY_s0; int Y0_x1 = Y0_s1 + SNY_s1; int Y0_q1 = Y0_x1 + SNY_x1; int Y0_r1 = Y0_q1 + SNY_q1; int SNY = Y0_r1 + SNY_r1; /* Create strips image and copy the raw strips: */ float_image_t *SIMG = float_image_new(NC, SNX, SNY); pst_kodak_copy_rectangle(IMG, P, xlo, xhi, ylo_r1, yhi_r1, 0, Y0_r1, SNX, SNY_r1, SIMG); pst_kodak_copy_rectangle(IMG, P, xlo, xhi, ylo_s1, yhi_s1, 0, Y0_s1, SNX, SNY_s1, SIMG); pst_kodak_copy_rectangle(IMG, P, xlo, xhi, ylo_s0, yhi_s0, 0, Y0_s0, SNX, SNY_s0, SIMG); pst_kodak_copy_rectangle(IMG, P, xlo, xhi, ylo_r0, yhi_r0, 0, Y0_r0, SNX, SNY_r0, SIMG); /* Compute the extrapolated and quotient strips: */ pst_kodak_copy_rectangle(IMG, P, xlo, xhi, ylo_r1, yhi_r1, 0, Y0_q1, SNX, SNY_q1, SIMG); pst_kodak_copy_rectangle(IMG, P, xlo, xhi, ylo_r0, yhi_r0, 0, Y0_q0, SNX, SNY_q0, SIMG); /* ---------------------------------------------------------------- */ void pst_kodak_add_global_term(pst_kodak_basis_t B, int n, double A[], double b[]); /* Adds to the {n × n} least-squares system {A,b} a term that forces the sample value 1.0 to be mapped to the linear intensity value 1.0. */ /* Add a term that forces 1.0 to map to 1.0: */ /* pst_kodak_add_global_term(B, n, A, b); */ void pst_kodak_add_global_term(pst_kodak_basis_t B, int n, double A[], double b[]) { /* Add to the quadratic functional the term {G = 0.5*h(1.0)^2} This means adding to each component {p} of the gradient {(A z - b)[p]} the term {(SUM{q: z[q]*bas[q](1)})*bas[p](1)} That is, we add {bas[q](1)*bas[p](1)} to each {A[p,q]}, and add nothing to {b[p]}. */ int p, q; /* Evaluate the basis functions at 1: */ double bas1[n]; for (p = 0; p < n; p++) { bas1[p] = pst_kodak_eval_basis(1.0, p, B); } for (p = 0; p < n; p++) { for (q = 0; q < n; q++) { A[p*n + q] += bas1[p]*bas1[q]; } } } /* ---------------------------------------------------------------- */ /* B-spline basis of degree 2, continuity 1: */ if (iv == k) { return 0.5*(1-r)*(1-r); } else if (iv == k-1) { return 0.5 + r*(1-r); } else if (iv == k-2) { return 0.5*r*r; } else { return 0.0; } /* ---------------------------------------------------------------- */ /* ---------------------------------------------------------------- */ SEE junk-minPower.c /* ---------------------------------------------------------------- */ /* ---------------------------------------------------------------- */ #define pst_lamp_spec_radiusOnly_HELP \ "-lamp \\\n" \ " [ " pst_lamp_spec_ambient_HELP " \\\n" \ " | " pst_lamp_spec_radius_HELP " \\\n" \ " | " pst_lamp_spec_wall_HELP " ]" \ #define pst_lamp_spec_radiusOnly_parameters_list_INFO \ "\"radius\", \"wall\", or \"ambient\"" #define pst_lamp_spec_radiusOnly_INFO \ "Begins the description of a light source. May be" \ " followed by " pst_lamp_spec_radiusOnly_parameters_list_INFO "." #define pst_lamp_spec_radiusOnly_HELP_INFO \ " -lamp {LAMP_PARAMETERS} \n" \ " " pst_lamp_spec_radiusOnly_INFO "\n" \ "\n" \ pst_lamp_spec_params_radiusOnly_HELP_INFO #define pst_lamp_spec_params_radiusOnly_HELP_INFO \ pst_lamp_spec_radius_HELP_INFO " Excludes the" \ " \"ambient\" and \"wall\" parameters.\n" \ "\n" \ pst_lamp_spec_wall_HELP_INFO " Excludes the" \ " \"radius\" and \"ambient\" parameters.\n" \ "\n" \ pst_lamp_spec_ambient_HELP_INFO " Excludes the" \ " \"radius\" and \"wall\" parameters." /* ---------------------------------------------------------------- */ void pst_flt_solve_lsq_system_nonneg(int NS, double Ab_s[], double u_s[]) { demand(FALSE, __FUNCTION__ " not implemented yet"); /* Subset of active variables: */ int NA = 0; /* Number of active variables. */ int six[NS]; /* Variable {u_a[ia]} is {u_s[six[ia]]} for {ia=0..NA-1}. */ int aix[NS]; /* Variable {u_s[is]} is {u_a[aix[is]]}, if {aix[is] >= 0}. */ /* Inactive variables are clamped at 0, active variables are non-negative. */ /* Start with all variables off, turn them on one at a time. */ { int is; for (is = 0; is < NS; is++) { u_s[is] = 0; } } /* Iterate until solution is satisfactory: */ while (TRUE) { /* In the current solution {u_s}, all elements are non-negative. */ /* Rebuild active set with all variables that are strictly positive: */ NA = 0; for (is = 0; is < NS; is++) { if (u_s[is] > 0) { /* Lamp is active: */ assert(NA < NS); six[NA] = is; aix[is] = NA; NA++; } } /* See if there is an inactive variable that should be activated: */ for (is = 0; is < NS; is++) { int ia = aix[is]; if (ia < 0) { /* Compute the derivative of the quadratic function at {u_s} {u_s[is]}: */ /* The gradient at a point {u} is {Au - b}: */ double dQdu = 0.0; for (js = 0; js < NS; js++) { dQdu = dQdu + Ab_s[is*(NS+1) + js]*u_s[js]; } dQdu = dQdu - Ab_s[is*(NS+1) + NS]; if (dQdu > 0.0) { /* Activate this lamp (it will improve the fit). */ assert(NA < NS); six[NA] = is; aix[is] = NA; NA++; /* Solve new system: */ } /* Check whether any active lamp has become negative: */ { int ia; for (ia = 0; ia < NA; ia++) { if (u_a[ia]) < 0.0 u_s[is] = (ia < 0 ? 0.0 : u_a[ia]); } } /* Build full solution vector {u_s} from active one {u_a}: */ { int is; for (is = 0; is < NS; is++) { int ia = aix[is]; u_s[is] = (ia < 0 ? 0.0 : u_a[ia]); } } pst_lamp_t *src = lmpv->el[is]; assert(src->pwr.nel == 1); if (u[ia] < 0) { /* Constrain lamp to zero: */ /* Set lamp's power to zero and forget about it: */ src->pwr.el[0] = 0.0; aix[is] = -1; } else { /* Set lamp's power to solution: */ src->pwr.el[0] = u[ia]; /* Keep lamp in new active subset: */ six[KA] = is; aix[is] = KA; KA++; } } int KA = 0; /* Number of variables retained in active subset. */ for (ia = 0; ia < NA; ia++) { int is = six[ia]; pst_lamp_t *src = lmpv->el[is]; assert(src->pwr.nel == 1); if (u[ia] < 0) { /* Constrain lamp to zero: */ /* Set lamp's power to zero and forget about it: */ src->pwr.el[0] = 0.0; aix[is] = -1; } else { /* Set lamp's power to solution: */ src->pwr.el[0] = u[ia]; /* Keep lamp in new active subset: */ six[KA] = is; aix[is] = KA; KA++; } } /* If no new constraints were introduced, we are done: */ if (KA == NA) { break; } /* Iterate with new active set: */ NA = KA; } /* Find best fit with currently active variables: */ if (NA > 0) { /* Extract from {Ab} the linear sub-system for the active variables only: */ double Ab_a[NA*(NA+1)]; /* Active parts of {A} and {b}, side-by-side. */ double u_a[NA]; /* Solution vector (powers of active variables). */ int ia; for (ia = 0; ia < NA; ia++) { int is = six[ia]; int ja; for (ja = 0; ja <= NA; ja++) { int js = six[ja]; Ab_a[ia*(NA+1)+ja] = Ab_s[is*(NS+1)+js]; } } /* Solve system: */ gsel_triangularize(NA, NA+1, Ab_a, TRUE); gsel_diagonalize(NA, NA+1, Ab_a, TRUE); gsel_normalize(NA, NA+1, Ab_a, TRUE); gsel_extract_solution(NA, NA+1, Ab_a, 1, u_a, TRUE); } } } } /* ---------------------------------------------------------------- */ /* Collect initial set of lamps to fit: */ /* ---------------------------------------------------------------- */ /* If the direction was undefined, supply an arbitrary default anyway: */ if ((dirAdjust) && (src->crad > -1.0) && (r3_L_inf_norm(&(src->dir)) == 0.0)) { fprintf(stderr, "%s: ", __FUNCTION__); fprintf(stderr, "undefined direction, setting it to (0,0,1)\n"); src->dir = (r3_t) {{ 0, 0, 1 }}; } /* ---------------------------------------------------------------- */ void pst_fit_light_simple_compute_lsq_mask ( float_image_t *IMG, /* Photo of object. */ float_image_t *NRM, /* Normal map of object. */ r3_t *dir, /* Estimated direction towards lamp. */ double crad, /* Co-sine of max angular radius of lamp. */ double minNormalZ, /* Ignore pixels where Z of normal is less than this. */ float_image_t *MSK /* (OUT) Weight mask for least-squares fitting. */ ); /* Fills the single-channel image {MSK} with a weight mask for least-squares fitting of the the simple light field model. Assumes {IMG} is the object's photo, {NRM} is its normal map, {dir} is the approximate direction of the lamp, and {crad} is the co-sine of the angular radius of a cone that contains the whole lamp. Excludes from {MSK} pixels where the normal's Z component is less than {minNormalZ}. Also excludes pixels where the sine of the normal and {dir} is greater than {crad} (unless {dir} is NULL or all zeros). The {MSK} image must be pre-allocated by the caller, with the same size as {IMG} and {NRM}. */ void pst_fit_light_simple_compute_lsq_mask ( float_image_t *IMG, float_image_t *NRM, r3_t *bdir, double crad, double minNormalZ, float_image_t *MSK ) { int NX = IMG->sz[1]; int NY = IMG->sz[2]; int x, y; demand((bdir == NULL) || (crad > 0.0), "lamp must be proper"); double srad = sqrt(fmax(0, 1 - crad*crad)); /* Sine of ang radius of source. */ for (y = 0; y < NY; y++) { for (x = 0; x < NX; x++) { r3_t nrm = pst_normal_map_get_pixel(NRM, x, y); double pix = float_image_get_sample(IMG, 0, x, y); double msk = 0.0; bool_t valid = (pix > 0.0) && (r3_L_inf_norm(&nrm)); if (valid) { bool_t shadowed = (bdir == NULL ? FALSE : r3_dot(&nrm, bdir) < srad); bool_t steep = (nrm.c[2] < minNormalZ); if ((! shadowed) && (! steep)) { msk = 1.0; } } float_image_set_sample(MSK, 0, x, y, msk); } } } /* ---------------------------------------------------------------- */ /* The procedure requires the scene's photo {IMG}, its normal map {NRM}, and a single-channel mask image {MSK} that specifies the pixels where the fitting is to be done, and their weights in the least-squares quadratic functional. */ /* ---------------------------------------------------------------- */ /* Build least squares system for main and ambient lamps within {MSK}: */ r4x4_t A; /* Scalar product matrix. */ r4_t b; /* Right-hand side vector. */ pst_fit_light_build_lsq_system(IMG, NRM, MSK, 4, bas, fun, &(A.c[0][0]), b.c); /* Solve system: */ r4_t u; /* Basis coefficient vector. */ r4x4_t I; r4x4_inv(&A, &I); r4x4_map_col(&I, &b, &u); /* Extract {bdir}, {bpwr}, and {apwr} from the coeff vector {u}: */ bdir->c[0] = u.c[0]; bdir->c[1] = u.c[1]; bdir->c[2] = u.c[2]; (*bpwr) = r3_dir(bdir, bdir); (*apwr) = u.c[4]; double bas(int i, r3_t *nrm) { return (i < 3 ? nrm->c[i] : 1.0); } double fun(r3_t *nrm, double img) { return img; } /* ---------------------------------------------------------------- */ #define pst_fit_light_small_dir_change 1.0e-4 /* Do not reactivate inactive lamps if change in {src.dir} was this small. */ /* Flag for active lamps (inactive ones must get {pwr=0}): */ bool_t active[NS]; /* {active[i]} is TRUE iff lamp {i} is . */ for (i = 0; i < NS; i++) { active[i] = TRUE; } /* If direction of {src} changed, re-activate all lamps: */ if (dir_change > pst_fit_light_small_dir_change) { for (i = 0; i < NS; i++) { active[i] = TRUE; } } /* Reset original power of all lamps: */ for (i = 0; i < NS; i++) { srcv.el[i]->pwr.el[0] = (active[i] ? ipwr[i] : 0.0); } /* ---------------------------------------------------------------- */ /* THE SIMPLE LIGHT FIELD MODEL */ #define pst_fit_light_simple_model_INFO \ "A /simple light model/ consists of two terms: an isotropic /ambient field/" \ " term, with arbitrary power, and a single /distant round lamp/ of" \ " bounded extent, with arbitrary direction, angular radius, and power." void pst_fit_light_simple_cast(pst_light_t *lht, pst_lamp_t **bnd, pst_lamp_t **amb); /* Makes sure that the light field {lht} contains exactly two distant round lamps, one with bounded extent ({crad > -1}) and one ambient ({crad == -1}). Adds those lamps to {lht} if necessary. Fails if {lht} has two or more lamps of either type. If it succeeds, returns in {*bnd} and {*amb} pointers to the bounded and ambient lamps in {lht}. */ void pst_fit_light_multi_cast(pst_light_t *lht); /* Makes sure that the light field {lht} has (1) at least one lamp, (2) at most one ambient lamp, and (3) at most three wall lamps. */ void pst_fit_light_multi_cast(pst_light_t *lht) { pst_lamp_vec_t *lmpv = &(lht->lmpv); int NS = lmpv->nel; /* There must be at least one lamp: */ if (lmpv->nel == 0) { /* Add an ambient lamp: */ pst_light_add_single(lmpv, &NS, -1.0); } pst_lamp_vec_trim(lmpv, NS); /* Count ambient lamps and wall lamps: */ int nwal = 0, namb = 0; int i; for (i = 0; i < NS; i++) { pst_lamp_t *src = lmpv->el[i]; if (src->crad <= -1.0) { namb++; } else if (src->crad == 0.0) { nwal++; } } demand(namb <= 1, "cannot fit two or more ambient lamps"); demand(nwal <= 3, "cannot fit four or more wall lamps"); } /* ---------------------------------------------------------------- */ /* Make a mask for the full image, minus parts that are too steep: */ pst_fit_light_simple_compute_lsq_mask(IMG, NRM, NULL, -1.0, minNormalZ, MSK); /* ---------------------------------------------------------------- */ /* pst_fit_light_single_lsq: */ /* The procedure requires a single-channel mask {MSK} which identifies those pixels where the fitting is to be done. The mask must be non-zero only in valid pixels (where {IMG} and {NRM} are non-zero), and only where the object is fully illuminated by the bounded lamp (i.e. excluding any parts that are in shadow, projected OR PROPER). The mask value {MSK[p]} is used as a weight in the quadratic functional used in the fitting. The mask should be guessed and/or adjusted iteratively by the caller. */ /* ---------------------------------------------------------------- */ /* pst_fit_light_single_iterative: */ /* If the parameter {MSK} is not NULL, the procedure stores in channel 0 of {*MSK} a mask that shows the weight of each pixel in the least-squares fitting step performed. The weight of each pixel lies between 0 and 1, and is nonzero only if the pixel was used for the fitting. The image {MSK} must be allocated by the caller, with at least one channel, and the same size as {IMG}. */ /* ---------------------------------------------------------------- */ #define pst_fit_light_dir_and_amb_together FALSE /* TRUE tells {pst_fit_light_single_iterative} to use a single least-squares step for all 4 parameters, at each iteration; FALSE says to use two least-square steps per iteration, one solving for the lamp's direction, the other for the powers (lamp and rest). */ if (pst_fit_light_dir_and_amb_together) { /* Single-step least squares (bounded and ambient lamps together): */ /* Make a mask for the parts of the gauge that are illuminated in full: */ pst_fit_light_simple_compute_lsq_mask(IMG, NRM, bdir, ctot, minNormalZ, MSK); /* Do least squares in that region for {bdir}, {bpwr}, {apwr}: */ pst_fit_light_single_lsq(IMG, NRM, MSK, bdir, bpwr, apwr); } else { /* Two-step least squares (first bounded lamp direction, then main/apwr powers). */ } /* ---------------------------------------------------------------- */ void pst_fit_light_single_direction_lsq ( float_image_t *IMG, float_image_t *NRM, pst_lamp_t *src, /* Lamp of {lht} to adjust. */ pst_light_t *lht, /* Light model. */ pst_lamp_t *src, /* Lamp of {lht} to adjust. */ float_image_t *MSK, /* Mask for least-squares error measure. */ float_image_t *MSK, r3_t *sdir, /* (OUT) Best-fitting direction for {src}. */ double *spwr /* (OUT) Best-fitting power for {src}. */ ); /* Returns in {*sdir} and {*spwr} a direction and power for {src} that would produce the best best fit between {IMG} and the white Lambertian shading model {F(lht,NRM)}, in the least-squares sense. Assumes that the remaining lamps are not to be adjusted. The procedure requires a single-channel mask {MSK} which identifies those pixels where the fitting is to be done. The mask must be non-zero only in valid pixels (where {IMG} and {NRM} are non-zero), and only where the object is fully illuminated by the bounded lamp (i.e. excluding any parts that are in shadow, projected OR PROPER). The mask value {MSK[p]} is used as a weight in the quadratic functional used in the fitting. The mask should be guessed and/or adjusted iteratively by the caller. */ #define pst_fit_light_single_direction_lsq_INFO \ " The procedure {pst_fit_light_single_direction_lsq}" \ " uses the following fact: at any pixel {p} where surface" \ " where surface is white Lambertian and the lamp is fully visible from" \ " it, the apparent lightness predicted by the simple light" \ " model is a linear function{F(NRM[p])} of the coordinates" \ " of the normal direction vector {NRM[p]}. The three unknown" \ " coefficients of the function {F} are the three coordinates" \ " of the lamp's direction times its power.\n" \ "\n" \ " Therefore, the procedure applies affine least squares fitting" \ " of {F(NRM[p])} to {IMG[p]} (over the fully" \ " lighted parts of the image, as indicated by the mask {MSK}); and" \ " then recovers the lamp's direction and power from the three" \ " coefficients of the fitted function {F}." #define pst_fit_light_single_intensities_by_least_squares_INFO \ " The procedure {pst_fit_light_simple_intensities_by_least_squares}" \ " uses the following fact: if the power {APWR} of the ambient lamp" \ " is known, then, at any pixel {p} where surface" \ " where surface is white Lambertian and the lamp is fully visible from" \ " it, the apparent lightness predicted by the simple light" \ " model is a linear function{F(NRM[p])} of the coordinates" \ " of the normal direction vector {NRM[p]}. The three unknown" \ " coefficients of the function {F}" \ " are the three coordinates of the lamp's direction times its" \ " power.\n" \ " Therefore, the procedure applies linear least squares fitting" \ " of {F(NRM[p])} to {IMG[p]-APWR} (over the fully" \ " lighted parts of the image, as indicated by the mask {MSK}); and" \ " then recovers the bounded lamp's direction and power, as" \ " well as the ambient lamp's power, from the four coefficients" \ " of the fitted function {F}." void pst_fit_light_simple_intensities_by_least_squares ( float_image_t *IMG, float_image_t *NRM, r3_t *bdir, double crad, float_image_t *MSK, double *bpwr, /* (OUT) Intensity of bounded lamp. */ double *apwr /* (OUT) Intensity of ambient lamp. */ ); /* Returns in {*bpwr} and {*apwr} the intensities of a bounded lamp and an ambient lamp that best fit {IMG} and {NRM}, assuming that the lamp has a known fixed direction {bdir} and cosine-radius {crad} The procedure requires a single-channel mask {MSK} which identifies those pixels where the fitting is to be done. The mask must be non-zero only in valid pixels (where {IMG} and {NRM} are non-zero). It should exclude projected shadows and highlights, but otherwise should include the whole surface, even parts in proper shadow on in the terminator zone. The mask value {MSK[p]} is used as a weight in the quadratic functional used in the fitting. */ /* ---------------------------------------------------------------- */ void pst_get_simple_model_params ( pst_light_t *lht, double *crad, r3_t **dir, double **pwr, double **amb ); /* Makes sure that {lht} is indeed a simple light model, consisting of a single bounded lamp {src} and a single ambient lamp {amb}, both monochromatic. Then sets {(crad} to {src.crad}, {*dir} to the address of {src.dir}, {*pwr} to the address of {src.pwr[0]}, and {*amb} to the address of {amb.pwr[0]}. In all three cases, the addresses point to the relevant fields inside the {lht} data structure. */ /* ---------------------------------------------------------------- */ /* fprintf(stderr, " max intensity = %8.5f\n", valMax); */ /* fprintf(stderr, " min intensity = %8.5f\n", valMin); */ /* fprintf(stderr, " max normal = "); */ /* r3_print(stderr, &nrmMax); */ /* fprintf(stderr, "\n"); */ /* fprintf(stderr, " coeffs = "); */ /* r4_print(stderr, &u); */ /* fprintf(stderr, "\n"); */ /* ---------------------------------------------------------------- */ /* COEFFICIENT VECTOR OF A SIMPLE LIGHT FIELD Actually, the fitting procedures encode the light field parameters as a single four-vector {u}, such that the first three components are the unit direction vector {dir} times the lamp's intensity {pwr}; while the fourth component is the ambient light intensity {amb}. */ void pst_simple_light_split_coeff_vector(r4_t *u, r3_t *dir, double *pwr, double *amb); /* Splits the coefficient {u} of a simple monochromatic light field model into the separate lamp parameters {src,pwr} and the ambient intensity {amb}. */ r4_t pst_simple_light_pack_coeff_vector(r3_t *dir, double pwr, double amb); /* Combines the parameters {dir,pwr} of single distant compact lamp and the ambient intensity {amb} into the coefficient vector {u} of a simple monochromatic light field model. */ void pst_simple_light_split_coeff_vector(r4_t *u, r3_t *dir, double *pwr, double *amb) { dir->c[0] = u->c[0]; dir->c[1] = u->c[1]; dir->c[2] = u->c[2]; (*pwr) = r3_dir(dir, dir); (*amb) = u->c[3]; } r4_t pst_simple_light_pack_coeff_vector(r3_t *dir, double pwr, double amb) { r4_t u; u.c[0] = pwr*dir->c[0]; u.c[1] = pwr*dir->c[1]; u.c[2] = pwr*dir->c[2]; u.c[3] = amb; return u; } /* ---------------------------------------------------------------- */ /* Compute projection matrices needed for normal computation: */ r3x3_t xym_to_uvm, uvw_to_xyz; pst_geom_sphere_view_matrices(geo, &xym_to_uvm, &uvw_to_xyz); fprintf(stderr, "entering %s\n", __FUNCTION__); fprintf(stderr, " source:\n"); pst_source_spec_write(stderr, src); fprintf(stderr, "\n"); fprintf(stderr, "entering %s\n", __FUNCTION__); for (i = 0; i < NL; i++) { pst_source_t *src = &(svec->el[i]); fprintf(stderr, " source %d:\n", i); pst_source_spec_write(stderr, src); fprintf(stderr, "\n"); } void pst_locate_sphere_compute_steps ( double ctrAdj, /* Max amount of adjustment in center coordinates. */ double radAdj, /* Max amount of adjustment in radius. */ double strAdj, /* Max amount of adjustment in stretch length. */ double adjStep, /* Ideal adjStep between trial parameter values. */ int maxTrials, /* Max number of trial parameter configs. */ int *KCP, /* (OUT) Number of trial values for center coords. */ int *KRP, /* (OUT) Number of trial values for radius. */ int *KEP /* (OUT) Number of trial values for stretch length. */ ); /* Computes the number of trial values for coordinates ({KC}), for the radius ({KR}), and for the stretch length ({KE}). These counts are always odd and positive. If any of {ctrAdj}, {radAdj}, and {strAdj} are positive, the cprresponding count is at least 3. Tries to use the specified {adjStep} in each parameter, but may use a larger step in order to keep the total number of trials below {maxTrials}, or at least close to it. */ void pst_locate_sphere_compute_steps ( double ctrAdj, /* Max amount of adjustment in center coordinates. */ double radAdj, /* Max amount of adjustment in radius. */ double strAdj, /* Max amount of adjustment in stretch length. */ double adjStep, /* Ideal step between trial parameter values. */ int maxTrials, /* Max number of trial parameter configs. */ int *KCP, /* (OUT) Number of trial values for center coords. */ int *KRP, /* (OUT) Number of trial values for radius. */ int *KSP /* (OUT) Number of trial values for stretch length. */ ) { /* Adjust {maxTrials} if less than the minimum needed for proper sampling: */ int minTrials = (ctrAdj > 0 ? 9 : 1)*(radAdj > 0 ? 3 : 1)*(strAdj > 0 ? 3 : 1); if (maxTrials < minTrials) { maxTrials = minTrials; } /* Compute number of values to try for each parameter. */ int KC = 1 + 2*(int)ceil(ctrAdj/adjStep); /* Num of trials in {ctr} X and Y. */ int KR = 1 + 2*(int)ceil(radAdj/adjStep); /* Num of trials in {rad}. */ int KS = 1 + 2*(int)ceil(strAdj/adjStep); /* Num of trials in {rad}. */ /* Try to honor the {maxTrials} limit: */ int KTot = KC*KC*KR*KS; if (KTot > maxTrials) { /* Adjust trial counts to ensure min 3 if varying, close to {maxTrials} total: */ int KVar = 2*(ctrAdj > 0.0) + (radAdj > 0.0) + (strAdj > 0.0); double scale = pow((double)maxTrials/(double)KTot, 1.0/(double)KVar); if (KC > 3) { KC = 1 + 2*(int)fmax(scale*(KC-1)/2 + 0.5, 1.0); } if (KR > 3) { KC = 1 + 2*(int)fmax(scale*(KR-1)/2 + 0.5, 1.0); } if (KS > 3) { KS = 1 + 2*(int)fmax(scale*(KS-1)/2 + 0.5, 1.0); } } /* Check and return results: */ assert((KC > 0) && (KC % 2 == 1)); (*KCP) = KC; assert((KR > 0) && (KR % 2 == 1)); (*KRP) = KR; assert((KS > 0) && (KS % 2 == 1)); (*KSP) = KS; } /* Compute actual steps to use in each parameter: */ double ctrStep = (KC == 1 ? 0.0 : 2*ctrAdj/(KC-1)); /* Actual step in {centerX,centerY} */ double radStep = (KR == 1 ? 0.0 : 2*radAdj/(KR-1)); /* Actual step in {radius} */ double strStep = (KS == 1 ? 0.0 : 2*strAdj/(KS-1)); /* Actual step in {strLen} */ /* ---------------------------------------------------------------------- */ /* pst_signature.c */ void func_peso(pixel* px, double_vec_t* peso, int canal); float dist(double_vec_t* u, double_vec_t* v); void func_peso(pixel* px, double_vec_t* peso ,int canal){ int i; float soma_cor = 0.0001; for(i = 0;i< peso->ne;i++){ float corAtual = ((float)px[i].canal[canal]); corAtual = corAtual*corAtual; soma_cor = soma_cor + corAtual; } soma_cor = sqrt(soma_cor); for(i = 0;i< peso->ne;i++){ float corAtual = ((float)px[i].canal[canal]); peso->e[i] = corAtual/soma_cor; } } /* ---------------------------------------------------------------------- *] /* pst_img_graph */ /* void pst_img_graph_remove_paralel_edges(pst_img_graph_t *g) */ /* { */ /* int32_t i; */ /* for (i = 0; i < g->m ; i++) */ /* { */ /* oct_arc_t e = g->edge[i]; */ /* if (e == oct_arc_NULL) continue; */ /* int32_t dst_e = pst_img_graph_get_edge_origin(g,oct_sym(e)); */ /* oct_arc_t a = oct_onext(e); */ /* int32_t count_repeated = 0; */ /* double sumW = pst_img_graph_get_edge_weight(g,e); */ /* double sumD = pst_img_graph_get_edge_delta(g,e)*sumW; */ /* while(a != e) */ /* { */ /* int32_t dst_a = pst_img_graph_get_edge_origin(g,oct_sym(a)); */ /* oct_arc_t a_next = oct_onext(a); */ /* if (dst_a == dst_e) */ /* { */ /* double wa = pst_img_graph_get_edge_weight(g,a); */ /* double da = pst_img_graph_get_edge_delta(g,a); */ /* sumW+=wa; */ /* sumD+=wa*da; */ /* count_repeated++; */ /* pst_img_graph_edge_remove(g,a); */ /* } */ /* a = a_next; */ /* } */ /* if (count_repeated > 0) */ /* { */ /* pst_img_graph_set_edge_delta(g,a,sumD); */ /* pst_img_graph_set_edge_weight(g,a,sumW); */ /* */ /* } */ /* } */ /* } */ /* ...................................................................... */ /* d = d/(wtot*wtot); */ /* double w; */ /* double we = w_edge[count]; */ /* double wf = w_edge[(count+1)%n_neighbours]; */ /* assert((we > 0) && (wf > 0)); */ /* double w_default = (we*wf)/(wtot); */ /* if (merge_diagonals) */ /* { */ /* double wd = w_edge[(count-1 + n_neighbours)%n_neighbours]; */ /* double wg = w_edge[(count+2)%n_neighbours]; */ /* w_default*= (1+ 0.5*(wg/wtot)*((we/wf) +1.0) + 0.5*(wd/wtot)*((wf/we)+1.0) ); */ /* } */ /* ...................................................................... */ /* int32_t *ix = (int32_t *)notnull(malloc((NXY_Z)*sizeof(int32_t)), "no mem"); */ /* { */ /* int32_t xy; */ /* for (xy = 0; xy < NXY_Z; xy++) */ /* { */ /* ix[xy] = -1; */ /* } */ /* int32_t count_idx = 0; */ /* for (xy = 0; xy n; xy++) */ /* { */ /* if ( ind_ix[xy] >= 0) */ /* { */ /* int32_t k = ind_ix[xy]; */ /* int32_t x,y; */ /* pst_img_graph_get_vertex_image_indices(&(g->vertex[xy]),NX_Z,NY_Z,&x,&y); */ /* */ /* col[count_idx] = x; */ /* row[count_idx] = y; */ /* ix[g->vertex[xy].id] =k; */ /* count_idx++; */ /* } */ /* } */ /* assert(count_idx == N); */ /* */ /* } */ /* ...................................................................... */ /* fprintf(stderr,"Checking consistency..."); */ /* for (i = 0; i < g->m; i++) */ { /* oct_arc_t e = g->edge[i].edge; */ /* if ( e != oct_arc_NULL) */ { /* int32_t ind = pst_img_graph_get_dir_edge_num(e); */ /* oct_arc_t e_eqv = ref_tab[ind]; */ /* */ /* assert(oct_onext(e_eqv) == linha(oct_onext(e)) ); */ /* */ /* assert(oct_onext(oct_sym(e_eqv)) == linha(oct_onext(oct_sym(e))) ); */ /* */ /* } */ /* } */ /* fprintf(stderr,"OK\n"); */ /* ---------------------------------------------------------------------- */ /* pst_virtual_gauge.{h,c} */ frgb_t pst_vitual_gauge_color_at_point ( r2_t *p, ellipse_crs_t *E, r3_t *anrm, frgb_t *albedo, int32_t NC, pst_virtual_gauge_shading_func_t* shading ); /* Returns the apparent color {val} of the gauge or background at point {p} of the synthetic image */ frgb_t pst_vitual_gauge_color_at_point ( r2_t *p, ellipse_crs_t *E, r3x3_t *VM, frgb_t *albedo, int32_t NC, pst_virtual_gauge_shading_func_t *shading ) { bool_t debug = FALSE; /* Subtract the center and undo the stretching: */ r2_t q = ellipse_crs_relative_coords(E, p); double r2 = r2_norm_sqr(&q); double val; if(r2 <= 1.0) { /* Compute the view-relative normal {vnrm} to the sphere: */ r3_t vnrm = (r3_t){{ q.c[0], q.c[1], sqrt(1 - r2) }}; /* View-relative normal */ /* Compute the absolute normal {anrm}: */ r3_t anrm; r3x3_map_row(&vnrm, VM, &anrm); val = shading(&anrm); if (debug) { fprintf(stderr," normal at q = ( %lf, %lf )\n", q.c[0], q.c[1]); r3_print(stderr, &anrm); fprintf(stderr," value = %lf\n", val); } } else { val = 0; } assert(isfinite(val)); frgb_t res; for (int32_t c = 0; c < NC; c++) { res.c[c] = (float)(val*albedo.c[c]); } return res; } /* ---------------------------------------------------------------------- */ /* pst_lamp.h */ /* If {*NC} is NULL, or the number of components {KC} in the parsed {COLOR} is exactly one, ignores {NC}. Otherwise, if {*NC} is negative, sets {*NC=KC} Otherwise demands that {KC} be equal to {*NC}. */ /* ---------------------------------------------------------------------- */ /* pst_light.{h,c} */ void pst_light_regularize_channels(pst_light_t *lht, int32_t NC); /* Makes sure that all lamps have {NC} channels in their power vectors. Vectors with zero elements are expanded {NC} elements and filled with {1/NS}, where {NS} is the number of lamps. Vectors with one channel are expanded to {NC} by replicating the first element. Vectors with {NC} elements are left untouched. A vector with any other length causes the procedure to bomb out. !!! Obsolete given that power is now an {frgb_t} !!! */ void pst_light_regularize_channels(pst_light_t *lht, int32_t NC) { pst_lamp_vec_t *lmpv = &(lht->lmpv); uint32_t NS = lmpv->ne; if (NS > 0) { double defpwr = 1.0/NS; for (uint32_t i = 0; i < NS; i++) { pst_lamp_t *src = lmpv->e[i]; pst_double_vec_regularize(&(src->pwr), NC, defpwr); /* The regularization must have succeded: */ assert(src->pwr.ne == NC); } } } /* ---------------------------------------------------------------------- */ /* pst_virtual_gauge.h */ float_image_t* pst_virtual_gauge_image_error ( float_image_t* img_in, float_image_t* img_plt, float_image_t* img_mk ); float_image_t* pst_shading_make_image_error ( float_image_t* img_in, float_image_t* img_plt, float_image_t* img_mk ) { int32_t NC = (int32_t)img_in->sz[0]; int32_t NX = (int32_t)img_in->sz[1]; int32_t NY = (int32_t)img_in->sz[2]; float_image_t* img_err = float_image_new(NC, NX, NY); for(int32_t c = 0; c < NC; c++){ for(int32_t x = 0; x < NX; x++){ for(int32_t y = 0; y < NY; y++){ double val1 = float_image_get_sample(img_in, c, x, y); double val2 = float_image_get_sample(img_plt, c, x, y); double mskval = float_image_get_sample(img_mk, 0, x, y); double mdif = (val1-val2)*mskval; float_image_set_sample(img_err, c, x, y, (float)mdif); } } } return img_err; }