/* See {multifok_sampling.h}. */ /* Last edited on 2024-12-09 06:25:07 by stolfi */ #include #include #include #include #include #include #include #include #include #include #include #include void multifok_sampling_choose_pixel_sampoints_and_weights ( uint32_t HS, uint32_t *NS_P, r2_t **uSmp_P, double **wSmp_P, bool_t verbose ); /* Generates an array {uSmp[0..NS-1]} of 2D vectors and an array {wSmp[0..NS-1]} of weights, where {NS = (2*HS+1)^2}. The vectors {uSmp[0..NS-1]} are a subset of the regular orthogonal grid of size {2*HS+1} by {2*HS+1}, symmetric about {(0,0)}, with spacing {1/(HS+1)} along each coordinate. Thus the max abs value of each coordinate is {HS/(HS+1)}. The weights are a 2D windowing function that is the product of a 1D Hann (raised cosine) window function on each coordinate. The vectors and weights form a partition or unit if replicated over the plane with stride 1.0 along each coordinate. The vectors and weights are sorted so that {|uSmp[k]|} is increasing with {k}. Thus {uSmp[0]} is always {(0,0)} and {wSmp[0]} is 1. Returns {NS,uSMP,wSmp} in {*NS_P,*uSmp_P,*wSmp_P}. */ void multifok_sampling_choose_ray_tilts_and_weights ( uint32_t KR_min, uint32_t NS, uint32_t *NR_P, r2_t **tRay_P, double **wRay_P, bool_t verbose ); /* Defines {NR} as {1} if {KR_min} is zero, otherwise as the smallest multiple of {NS} that is at least {KR_min*NS} and is the number of points in a digital disk. Alocates and fills {uRay[0..NS-1]}, {wRay[0..NS-1]}. Returns them in {*NS_P}, {*uRay_P}, {*wRay_P}. */ void multifok_sampling_get_grid_points_in_quadrant ( uint32_t HP, i2_t **pq_P, uint32_t *NPQ_P, bool_t verbose ); /* Allocates and fills a list {pq[0..NPQ-1]} of the {NPQ=HP*(HP+1)} integer grid points in the first quadrant with coordinates up to {HP}; specifically, {(ix,iy)} with {ix} in {1..HP} and {iy} in {0..HP}, sorted by distance from the origin. */ void multifok_sampling_get_digital_disk_points_in_quadrant ( uint32_t NS, uint32_t KP_min, i2_t **pq_P, uint32_t *NPQ_P, bool_t verbose ); /* Allocates and fills a list {pq[0..NPQ-1]} with the {NPQ} integer grid points in the first quadrant (namely points {(ix,iy)} with {ix >= 1) and {iy >= 0}) which are closest to the origin, sorted by distance from the origin. The number {NPQ} will be such that {NP=1+4*NPQ} is a multiple of {NS}, and {NP/NS} is at least {KP_min}. Both {KP_min} and {NS} musy be positive, and {NS} must be odd. */ void multifok_sampling_sort_points_by_norm(uint32_t NP, i2_t p[]); /* Sorts the list of points {p[0..NP-1]} by increasing distance from origin. Specifically, according to {multifok_sampling_compare_points}. */ void multifok_sampling_check_point_order(uint32_t NP, i2_t p[]); /* Checks whether the points {p[0..NP-1]} are sorted in incrasing order of distance from the origin. */ int32_t multifok_sampling_compare_points(const void *a, const void *b); /* Suitable for sorting routines like {qsort}. Assumes {a} and {b} point to integer pair ({i2_t}) records {pa} and {pb}, and returns {-1} or {+1} if {pa} should come before or after {pb}, respectively; or 0 if they are the same pair. Compares {pa} with {pb} by distance from origin. Breaks ties by {min(|x|,|y|)}, then by {|x|+|y|}, then by lex order. */ double *multifok_sampling_hann_weight_table(uint32_t HW, bool_t verbose); /* Allocates and fills a 1-dim table {w[0..NW-1]} of Hann (raised cosine) weights, where {NW = 2*HW + 1}. Element {w[HW]} will be 1.0, and {w[HW-i]} will be equal to {w[HW+i]}. The table will have the partition-of-unit property with stride {HW+1}; that is, {w[HW+1+i]+w[i]=1.0} for {i} in {0..HW-1}. */ void multifok_sampling_print_points_and_weights ( uint32_t NP, r2_t p[], char *pName, double w[], char *wName ); /* Prints {p0..NP-1]} and {w[0..NP-1]} to {stderr}. */ /* IMPLEMENTATIONS */ multifok_sampling_t *multifok_sampling_choose(uint32_t HS, uint32_t KR_min, bool_t verbose) { multifok_sampling_t *samp = talloc(1, multifok_sampling_t); multifok_sampling_choose_pixel_sampoints_and_weights (HS, &(samp->NS), &(samp->uSmp), &(samp->wSmp), verbose); multifok_sampling_choose_ray_tilts_and_weights (KR_min, samp->NS, &(samp->NR), &(samp->tRay), &(samp->wRay), verbose); return samp; } void multifok_sampling_free(multifok_sampling_t *samp) { free(samp->uSmp); free(samp->wSmp); free(samp->tRay); free(samp->wRay); free(samp); } void multifok_sampling_choose_pixel_sampoints_and_weights ( uint32_t HS, uint32_t *NS_P, r2_t **uSmp_P, double **wSmp_P, bool_t verbose ) { /* This implementation chooses the subsampling points as a regular orthogonal grid of {2*HS+1} by {2*HS+1}, and 2D Hann (raised cosine) weights. This ensures the partition of unity property for the sampling points. */ demand(HS >= 0, "invalid {HS}"); bool_t debug = TRUE; if (debug) { verbose = TRUE; } /* Allocate and fill the temporary arrays: */ uint32_t NW = 2*HS + 1; uint32_t NS = NW*NW; if (debug) { fprintf(stderr, "generating %d × %d = %d sampling vectors total (HS = %d)\n", NW, NW, NS, HS); } r2_t *uSmp = talloc(NS, r2_t); double *wSmp = talloc(NS, double); /* Define element 0: */ uint32_t ko = 0; /* Elements already filled: */ uSmp[ko] = (r2_t){{ 0,0 }}; wSmp[ko] = 1.0; ko++; if (HS > 0) { double *wHan = multifok_sampling_hann_weight_table(HS, verbose); uint32_t NSQ; i2_t *pq; multifok_sampling_get_grid_points_in_quadrant(HS, &pq, &NSQ, verbose); assert(NSQ > 0); assert(pq != NULL); assert(NS == 1 + 4*NSQ); /* Replicate the quadrant points into {uSmp[1..NS_max-1],wSmp[1..NS_max-1]: */ double scale = 1.0/(HS+1); for (uint32_t j = 0; j < NSQ; j++) { int32_t xj = pq[j].c[0]; assert(abs(xj) <= HS); int32_t yj = pq[j].c[1]; assert(abs(yj) <= HS); double wj = wHan[xj+(int32_t)HS]*wHan[yj+(int32_t)HS]; assert (ko + 4 <= NS); wSmp[ko] = wj; uSmp[ko] = (r2_t){{ +xj*scale, +yj*scale }}; ko++; wSmp[ko] = wj; uSmp[ko] = (r2_t){{ -yj*scale, +xj*scale }}; ko++; wSmp[ko] = wj; uSmp[ko] = (r2_t){{ -xj*scale, -yj*scale }}; ko++; wSmp[ko] = wj; uSmp[ko] = (r2_t){{ +yj*scale, -xj*scale }}; ko++; } free(pq); free(wHan); } assert(ko == NS); if (verbose) { fprintf(stderr, " generated %d sampling points and weights:\n", NS); multifok_sampling_print_points_and_weights(NS, uSmp, "uSmp", wSmp, "wSmp"); } (*NS_P) = NS; (*uSmp_P) = uSmp; (*wSmp_P) = wSmp; } void multifok_sampling_get_grid_points_in_quadrant ( uint32_t HP, i2_t **pq_P, uint32_t *NPQ_P, bool_t verbose ) { demand(HP >= 0, "invalid {HP}"); bool_t debug = TRUE; /* Allocate and fill the temporary arrays: */ uint32_t NW = 2*HP + 1; /* Generate the 2D samples and weights for the first quadrant: */ uint32_t NPQ = (NW*NW - 1)/4; if (debug) { fprintf(stderr, " generating %d grid points in first quadrant\n", NPQ); } i2_t *pq = talloc(NPQ, i2_t); uint32_t kq = 0; /* First-quadrant points generated so far. */ for (int32_t ix = 1; ix <= +HP; ix++) { for (int32_t iy = 0; iy <= +HP; iy++) { assert(kq < NPQ); pq[kq] = (i2_t){{ ix, iy }}; kq++; } } assert(kq == NPQ); multifok_sampling_sort_points_by_norm(NPQ, pq); (*pq_P) = pq; (*NPQ_P) = NPQ; } void multifok_sampling_choose_ray_tilts_and_weights ( uint32_t KR_min, uint32_t NS, uint32_t *NR_P, r2_t **tRay_P, double **wRay_P, bool_t verbose ) { uint32_t NRQ = 0; i2_t *pq = NULL; if (KR_min == 0) { if (verbose) { fprintf(stderr, " using one vertical ray for all sampling points\n"); } NRQ = 0; } else { demand(NS >= 1, "{NS} must be at least 1"); demand((NS % 2) == 1, "{NS} must be odd"); if (verbose) { fprintf(stderr, " trying to get at least %d rays for each of %d sampling points\n", KR_min, NS); } multifok_sampling_get_digital_disk_points_in_quadrant(NS, KR_min, &pq, &NRQ, verbose); } uint32_t NR = 1 + 4*NRQ; r2_t *tRay = talloc(NR, r2_t); double *wRay = talloc(NR, double); /* Define the first ray as straight down: */ uint32_t ko = 0; /* Count of total rays saved. */ tRay[ko] = (r2_t){{ 0, 0 }}; wRay[ko] = 1.0; ko++; if (NRQ > 0) { assert(NR >= NS); assert(NR % NS == 0); assert(NR/NS >= KR_min); assert(pq != NULL); /* Jitter the the quadrant points {tq[0..NRQ-1]} points then replicate to the other quadrants, appending them to {tRay[0..ko-1]} and computer their weights {wRay[0..ko-1]}. The jitter is applied before replication to keep the tilt pattern symmetrical about the origin. */ double jitter = 0.2; /* Max amount of jitter to apply. Should be less than 0.5. */ auto double jit(uint32_t k); /* Maps{k} to a pseudorandom function in {[-jitter _ +jitter]}, with roughly uniform distribution. */ for (uint32_t kk = 0; kk < 20; kk++) { fprintf(stderr, " jit = %14.10f\n", jit(kk)); } double rmax = sqrt((double)i2_norm_sqr(&(pq[NRQ-1]))); /* Max norm of accepted points. */ double sigma = rmax/3.0; /* Width of Gaussian bell. */ for (uint32_t iq = 0; iq < NRQ; iq++) { double xi = (double)pq[iq].c[0] + jit(iq); double yi = (double)pq[iq].c[1] + jit(iq); double ri = hypot(xi, yi); /* Jittered. */ double zi = ri/sigma; double wi = exp(-zi*zi/2); assert (ko + 4 <= NR); wRay[ko] = wi; tRay[ko] = (r2_t){{ +xi, +yi }}; ko++; wRay[ko] = wi; tRay[ko] = (r2_t){{ -yi, +xi }}; ko++; wRay[ko] = wi; tRay[ko] = (r2_t){{ -xi, -yi }}; ko++; wRay[ko] = wi; tRay[ko] = (r2_t){{ +yi, -xi }}; ko++; } assert(ko == NR); free(pq); /* Normalize the tilts to unit RMS radius: */ double sum_w = 0; double sum_w_r2 = 0; for (uint32_t ir = 0; ir < NR; ir++) { sum_w += wRay[ir]; double r2 = r2_norm_sqr(&(tRay[ir])); sum_w_r2 += wRay[ir]*r2; } assert(sum_w >= 1.0); double r_avg = sqrt(sum_w_r2/sum_w); assert(r_avg > 0.0); for (uint32_t ir = 0; ir < NR; ir++) { tRay[ir].c[0] /= r_avg; tRay[ir].c[1] /= r_avg; } double jit(uint32_t k) { double ran1 = sqrt(k + M_PI); ran1 = ran1 - floor(ran1); double ran2 = ran1*M_PI; ran2 = ran1 - floor(ran2); return jitter*ran2; } } if (verbose) { fprintf(stderr, " generated %d ray tilts and weights:\n", NR); multifok_sampling_print_points_and_weights(NR, tRay, "tRay", wRay, "wRay"); } /* Return results: */ (*NR_P) = NR; (*tRay_P) = tRay; (*wRay_P) = wRay; } void multifok_sampling_get_digital_disk_points_in_quadrant ( uint32_t NS, uint32_t KP_min, i2_t **pq_P, uint32_t *NPQ_P, bool_t verbose ) { bool_t debug = TRUE; demand(KP_min > 0, "invalid {KP_min}"); demand(NS >= 1, "invalid {NS}"); demand((NS % 2) == 1, "{NS} must be odd"); uint32_t NP_max = (KP_min == 0 ? 1 : 20001); /* A stupid upper bound for {NP}. */ uint32_t NPQ_max = 2*(NP_max - 1)/4; /* Upper bound for candidate points in quadrant. */ /* Generate rays in the first quadrant of the disk: */ if (debug) { fprintf(stderr, " generating up to %d candidate rays in first quadrant\n", NPQ_max); } i2_t *pq = talloc(NPQ_max, i2_t); /* Enumerate candidate integer points of quadrant in increasing order of {max(ix,iy)}: */ uint32_t mq = 0; /* First-quadrant tilt candidates generated so far. */ uint32_t kq = 0; /* First-quadrant tilts accepted so far. */ uint32_t HRQ = 0; /* {X} and {Y} of last generated candidate layer in quadrant. */ auto void get_next_cand(void); /* Assumes that the points {pq[0..mq]} are all the integer grid points in {1..HRQ} × {0..HRQ}. Assumes that the "accepted" points {pq[0..kq-1]} are not father from the origin than any of the points {pq[kq..mq-1]}. Appends to {pq[0..mq-1]} zero or more layers of candidate grid points starting at {pq[mq]}, and increments {mq}. Each layer consists of pairs {(1..xym,xym)} and {(xym,0..xym-1)}, starting with {xym=HRQ}. Increments {HRQ} as it adds each layer. Keeps adding layers until {HRQ} is definitely more than the distance from the origin to the last accepted point, plus one. Then sorts the candidates in the queue ({pq[kq..mq-1]}) by increasing distance from origin. This ensure that the next candidate in the queue {pq[kq]} exists and can be added to the "accepted" set without violating the above invariant. */ uint32_t NP; /* Total points in whole disk so far. */ while (TRUE) { NP = 1 + 4*kq; /* Total rays if we use {kq} per quadrant. */ /* Did we get the proper number of rays? */ if ((NP >= NS) && ((NP % NS) == 0) && (NP/NS >= KP_min)) { break; } /* Did we get too many rays already? */ demand(NP < NP_max, "failed to deermine a suitable set of rays"); /* We need to accept at least one more tilt per quadrant: */ if (debug) { uint32_t KP_old = NP/NS; uint32_t KP_new = (NP + 4)/NS; if ((KP_new >= KP_min) && (KP_new != KP_old)) { fprintf(stderr, " trying to use at least %d rays per sampling point ...\n", KP_new); } } get_next_cand(); assert(mq == HRQ*(HRQ+1)); /* Check order of points including last accepted one: */ uint32_t iniq = (kq == 0 ? 0 : kq-1); multifok_sampling_check_point_order(mq-iniq, &(pq[iniq])); /* Now accept one more point: */ kq++; } uint32_t NPQ = kq; /* Number of sample points actually retained in the first quadrant. */ if (debug) { fprintf(stderr, " retained %d candidate ray tilts in first quadrant\n", NPQ); } /* Seems to have converged: */ assert((NP % NS) == 0); uint32_t KP = NP/NS; assert(KP >= KP_min); (*pq_P) = pq; (*NPQ_P) = NPQ; return; void get_next_cand(void) { /* Get the radius of the last accepted point: */ double rprev = (kq == 0 ? 0.0 : sqrt((double)i2_norm_sqr(&(pq[kq-1])))); /* Make sure that we generated the required layer: */ while ((double)HRQ < rprev + 1.1) { /* Generate one more candidate layer: */ HRQ++; for (int32_t ixy = 0; ixy < HRQ; ixy++) { assert(mq + 2 <= NPQ_max); pq[mq] = (i2_t){{ ixy + 1, (int32_t)HRQ }}; mq++; pq[mq] = (i2_t){{ (int32_t)HRQ, ixy }}; mq++; } } assert(mq > kq); /* Re-sort {pq[kq..mq-1]} by distance: */ multifok_sampling_sort_points_by_norm(mq-kq, &(pq[kq])); } } void multifok_sampling_sort_points_by_norm(uint32_t NP, i2_t p[]) { qsort(p, NP, sizeof(i2_t), &multifok_sampling_compare_points); } void multifok_sampling_check_point_order(uint32_t NP, i2_t p[]) { bool_t debug = FALSE; for (uint32_t j = 0; j < NP; j++) { if (debug) { double rj = sqrt((double)i2_norm_sqr(&(p[j]))); fprintf(stderr, " p[%4d] = ", j); i2_gen_print(stderr, &(p[j]), "%+8d", " ( ", " ", " )"); fprintf(stderr, " norm = %15.9f\n", rj); } if (j > 0) { assert(multifok_sampling_compare_points(&(p[j-1]), &(p[j])) == -1); } } } int32_t multifok_sampling_compare_points(const void *a, const void *b) { i2_t *pa = ((i2_t *)a); uint64_t r2a = i2_norm_sqr(pa); i2_t *pb = ((i2_t *)b); uint64_t r2b = i2_norm_sqr(pb); if (r2a < r2b) { return -1; } else if (r2a > r2b) { return +1; } else { int32_t xa = pa->c[0], ya = pa->c[1]; int32_t xb = pb->c[0], yb = pb->c[1]; /* Try to break the tie by {min(|x|,|y|)} */ int32_t ma = (int32_t)imin(abs(xa),abs(ya)); int32_t mb = (int32_t)imin(abs(xb),abs(yb)); if (ma < mb) { return -1; } else if (ma > mb) { return +1; } else { /* Try to break the tie by {|x|+|y|} */ int32_t sa = abs(xa) + abs(ya); int32_t sb = abs(xb) + abs(yb); if (sa < sb) { return -1; } else if (sa > sb) { return +1; } else { /* Try to break the tie by lex order: */ if (xa < xb) { return -1; } else if (xa > xb) { return +1; } else if (ya < yb) { return -1; } else if (ya > yb) { return +1; } else { /* The pairs are equal: */ return 0; } } } } } void multifok_sampling_print_points_and_weights ( uint32_t NP, r2_t p[], char *pName, double w[], char *wName ) { for (uint32_t k = 0; k < NP; k++) { fprintf(stderr, " %4d %s = ", k, pName); r2_gen_print(stderr, &(p[k]), "%+9.6f", "( ", " ", " )"); if (w != NULL) { fprintf(stderr, " %s = %12.10f", wName, w[k]); } fprintf(stderr, "\n"); } } double* multifok_sampling_hann_weight_table(uint32_t HW, bool_t verbose) { uint32_t NW = 2*HW + 1; double *wHan = talloc(NW, double); wt_table_hann_fill(NW, 0.0, wHan, NULL); if (verbose) { for (uint32_t j = 0; j < NW; j++) { fprintf(stderr, " wHan[%d] = %10.8f\n", j, wHan[j]); } } assert(wHan[HW] == 1.0); /* Paranoia check of partition-of-unity property: */ for (uint32_t j = 0; j < HW; j++) { assert(fabs(wHan[j] + wHan[j+HW+1] - 1.0) < 1.0e-13); } return wHan; }