/* Last edited on 2024-12-26 15:01:38 by stolfi */ /* ------------------------------------------------------------------------ */ /* float_image_geostereo */ /* Smoothing weights ({WT0} = central pixel, {WT1} = adjacent): */ #define WT0 (0.425) #define WT1 (0.075) /* #define WTERR (2.0*sqrt(WT1*WT1*WT1*WT1 + 2.0*WT0*WT0*WT1*WT1 + WT0*WT0*WT0*WT0)) */ /* This may be a start for an optimization for a special case of {float_image_interpolate_grid_pixels} with {dx=dy=1.0} and {hy} integer. */ void float_image_geostereo_get_samples ( float_image_t *f, /* Pixel row buffer for image 1. */ uint32_t x, /* Reference column index. */ uint32_t y, /* Reference row index. */ double d, /* Horizontal displacement (pixels, fractional). */ uint32_t nwx, /* Window width. */ uint32_t nwy, /* Window height. */ float smp[] /* (OUT) Window sample buffer. */ ); /* Extracts samples from image {f} contained in a window centered on column {x+d} and row {y}. The sampling window will have {nwx} columns and {nwy} rows, centered at the specified indices. The window dimensions {nwx} and {nwy} must be odd. The displacement {d} can be fractional, in which case the window samples are obtained by smooth interpolation of the image samples. The output array {smp} must have enough space for all window samples; that is, {nwx*nw*NC} elements where {NC} is the number of channels in the input image. The sample on column {ix}, row {iy}, and channel {ic} of the window is stored in {smp[ic + NC*(ix + nwx*iy)]}. */ void float_image_geostereo_get_samples ( float_image_t *f, /* Pixel row buffer for image 1. */ uint32_t x, /* Reference column index. */ uint32_t y, /* Reference row index. */ double d, /* Horizontal displacement (pixels, fractional). */ uint32_t nwx, /* Window width. */ uint32_t nwy, /* Window height. */ float smp[] /* (OUT) Window sample buffer. */ ) { int32_t NC = (int32_t)f->sz[0]; int32_t NX = (int32_t)f->sz[1]; int32_t NY = (int32_t)f->sz[2]; int32_t rx = nwx/2; int32_t ry = nwy/2; /* Split {xd = x+d} into integer {xdi} and fraction {xdf}: */ double xd = (double)x + d; int32_t xdi = (int32_t)floor(xd); double xdf = xd - (double)xdi; /* Disregard very small fractions, maybe due to rounding: */ if (fabs(xdf) < 1.0e-14) { xdf = 0.0; } if (fabs(1.0 - xdf) < 1.0e-14) { xdi += 1; xdf = 0.0; } /* Check if there is any window pixel center inside the domain at all: */ bool_t inside_x = (xd+rx >= 0.0) && (xd-rx < (double)NX)); bool_t inside_y = (y+ry >= 0) && (y-ry < NY); if (! (inside_x && inside_y)) { /* Window is entirely outside the image: */ int32_t nsmp_tot = nwx*nwy*NC; for (int32_t k = 0; k < nsmp_tot; k++) { smp[k] = NAN; } } else if (xdf = 0.0) { /* Displacement is integer, just extract the samples: */ float_image_get_window_pixels(f, xdi, y, nwx, nwy, FALSE, smp); } else { /* Displacement is fractional, must interpolate with 4-tap cubic: */ int nsmp_row = NC*nwx; /* Samples to return per row. */ /* Find column range {xmin_raw..xmax_raw} that we need to get from {f}: */ int32_t xmin_raw = xdi - rx - 1; /* Column of first sample. */ int32_t xmax_raw = xdi + rx + 2; /* Column of last sample. */ int32_t nwx_raw = xmax_raw - xmin_raw + 1; /* Number of pixels to fetch from {f}: */ assert(nwx_raw == nwx+3); /* Paranoia. */ float smp_raw[nwx_raw]; /* Raw samples to be interpolated. */ int krow = 0; /* Index of first sample of next pixel row in {smp}. */ for(int32_t iy = -ry; iy <= +ry; iy++) { int yp = y + iy; if ((yp < 0) || (yp >= NY)) { /* This row is all {NAN}: */ int32_t ksmp = krow; for (int32_t icx = 0; icx < nsmp_row; icx++) { smp[ksmp] = NAN; ksmp++; } } else { for (int32_t c = 0; c < NC; c++) { /* Extract one row of samples, with {NAN} for out-of-bounds: */ float_image_get_sample_row(f, c, xmin_raw, xmax_raw, yp, smp_raw); /* Compute the desired samples: */ int32_t ksmp = krow + c; /* index of next sample of this channel in {smp}. */ for (int32_t sx = 0; sx < nwx; sx++) { float sa = smp_raw[ix + 0]; float sb = smp_raw[ix + 1]; float sc = smp_raw[ix + 2]; float sd = smp_raw[ix + 3]; smp[ksmp] = float_image_geostereo_interpolate(sa, sb, sc, sd, xdf); ksmp += NC; } } } krow += nsmp_row; } } } float float_image_geostereo_interpolate(float sa, float sb, float sc, float sd, int32_t rd) { if (rd == 1) { return (float)(- 0.1334*sa + 0.8169*sb + 0.3902*sc - 0.0507*sd); } else if (rd == 2) { return (float)(- 0.0507*sa + 0.3902*sb + 0.8169*sc - 0.1334*sd); } assert(0); } /* ----------------------------------------------------------------- */ int32_t rd = (d + 3000000) % 3; /* Fractional displacement in 1/3 pixs. */ int32_t id = (d - rd)/3; /* Displacement in whole pixels. */ int32_t nw = 0, ic, iy, ix; /* Index range for interpolation (rel to {id}): */ int32_t ixlo = -1; int32_t ixhi = (rd == 0 ? +1 : +2); for (iy = y-ry; iy <= y+ry; iy++) { bool_t iyok = ((iy >= 0) || (iy < NY)); for (ix = x+id-rx; ix <= x+id+rx; ix++) { bool_t ixok = ((ix+ixlo >= 0) || (ix+ixhi < NX)); for (ic = 0; ic < NC; ic++) { if (! (iyok && ixok)) { /* Out of bounds (partially or totally): */ w[nw] = NAN; } else if (rd == 0) { /* Exact index: */ w[nw] = float_image_get_sample(f, ic, ix, iy); } else { /* Fractional index, interpolate: */ float sa = float_image_get_sample(f, ic, ix - 1, iy); float sb = float_image_get_sample(f, ic, ix + 0, iy); float sc = float_image_get_sample(f, ic, ix + 1, iy); float sd = float_image_get_sample(f, ic, ix + 2, iy); w[nw] = float_image_geostereo_interpolate(sa, sb, sc, sd, rd); } nw++; } } /* ------------------------------------------------------------------------ */ /* float_image_read_pnm */ #define VIEW_GAMMA 2.200 #define VIEW_BIAS 0.010 /* Gamma and bias for viewing images on IBM-PC-like platforms. */ /* ------------------------------------------------------------------------ */ /* jspng_image_read */ /* Get automatic conversion of palette to RGB: */ if (color_type == PNG_COLOR_TYPE_PALETTE) { if (verbose) { fprintf(stderr, "requesting conversion from PALETTE to RGB\n"); } png_set_palette_to_rgb(pr); uint8_t old_color_type = color_type; color_type = png_get_color_type(pr, pi); if (color_type == PNG_COLOR_TYPE_PALETTE) { fprintf(stderr, "** {color_type} did not change from PALETTE\n"); } else if (color_type != PNG_COLOR_TYPE_RGB) { fprintf(stderr, "** {color_type} did not become RBG\n"); } assert(color_type == PNG_COLOR_TYPE_RGB); /* ??? I suppose ??? */ uint8_t old_bits = bits; bits = png_get_bit_depth(pr, pi); if (bits != old_bits) { fprintf(stderr, "!! {bits} changed from %u to %u\n", (uint32_t)old_bits, (uint32_t)bits); } } /* Extract sub-byte samples into separate bytes: */ if ((color_type == PNG_COLOR_TYPE_GRAY) && (bits < 8)) { if (verbose) { fprintf(stderr, "requesting expansion from %u bits per prixel to 8\n", (uint32_t)bits); } png_set_gray_1_2_4_to_8(pr); uint8_t old_bits = bits; bits = png_get_bit_depth(pr, pi); if (bits == old_bits) { fprintf(stderr, "* {bits} did not change\n"); } else if (bits != 8) { fprintf(stderr, "* {bits} became %d instead of 8\n", (uint32_t)bits); } } if (png_get_valid(pr, pi, PNG_INFO_tRNS)) { if (verbose) { fprintf(stderr, "adding alpha channel from \"tRNS\" chunk\n"); } if ((color_type == PNG_COLOR_TYPE_GRAY_ALPHA) || (color_type == PNG_COLOR_TYPE_RGB_ALPHA)) { fprintf(stderr, "** {color_type} already specifies ALPHA\n"); } png_set_tRNS_to_alpha(pr); uint8_t old_color_type = color_type; color_type = png_get_color_type(pr, pi); if (color_type == old_color_type) { fprintf(stderr, "** {color_type} did not change\n"); } else if ((color_type != PNG_COLOR_TYPE_GRAY_ALPHA) && (color_type != PNG_COLOR_TYPE_RGB_ALPHA)) { fprintf(stderr, "** {color_type} did not get ALPHA\n"); } assert((color_type == PNG_COLOR_TYPE_GRAY_ALPHA) || (color_type == PNG_COLOR_TYPE_RGB_ALPHA)); uint8_t old_bits = bits; bits = png_get_bit_depth(pr, pi); if (bits != old_bits) { fprintf(stderr, "!! {bits} changed from %u to %u\n", (uint32_t)old_bits, (uint32_t)bits); } } /* PNG files pack pixels of bit depths 1, 2, and 4 into bytes as small as they can, resulting in, for example, 8 pixels per byte for 1 bit files. This code expands to 1 pixel per byte without changing the values of the pixels:*/ if (bits < 8) { png_set_packing(pr); } /* PNG files have possible bit depths of 1, 2, 4, 8, and 16. All pixels stored in a PNG image have been "scaled" or "shifted" up to the next higher possible bit depth (e.g. from 5 bits/sample in the range [0,31] to 8 bits/sample in the range [0, 255]). However, it is also possible to convert the PNG pixel data back to the original bit depth of the image. This call reduces the pixels back down to the original bit depth: */ png_color_8p sBIT; if (png_get_sBIT(pr, pi, &sBIT)) { png_set_shift(pr, sBIT); } /* PNG files store 3-color pixels in red, green, blue order. This code changes the storage of the pixels to blue, green, red: */ if (color_type == PNG_COLOR_TYPE_RGB || color_type == PNG_COLOR_TYPE_RGB_ALPHA) { png_set_bgr(pr); } /* PNG files store RGB pixels packed into 3 or 6 bytes. This code expands them into 4 or 8 bytes for windowing systems that need them in this format: */ if (color_type == PNG_COLOR_TYPE_RGB) { png_set_filler(pr, filler, PNG_FILLER_BEFORE); } int32_t transforms = 0 | PNG_TRANSFORM_PACKING /* Expand 1, 2 and 4-bit samples to bytes. */ | PNG_TRANSFORM_EXPAND /* Expand palette, add alpha from tRNS chink. */ | PNG_TRANSFORM_SHIFT; /* Normalize pixels to the sBIT depth. */ png_read_png(pr, pi, NULL); /* ------------------------------------------------------------------------ */ /* float_image_transform.c */ void fitr_debug_direction(char *label, r2_t *v, char *tail); /* Prints the vector {v} to {stderr}. */ void fitr_debug_direction(char *label, r2_t *v, char *tail) { fprintf(stderr, " %s", label); fprintf(stderr, " = (%+8.5f %+8.5f)", v->c[0], v->c[1]); fprintf(stderr, "%s", tail); } /* Allocate scanlines of PNG image: */ png_bytep *png_arrayP = png_get_rows(pr, pi); /* -- INPUT TRANSFORMATION REQUESTS -------------------------------- */ /* Request expansion of 1,2,4-bit samples to one byte: */ /* This will leave each sample in the LOWER {bits} bits of each byte. */ if (bits < 8) { png_set_packing(pr); } /* Request conversion of palette colors to RGB or RGBA: */ if (color_type == PNG_COLOR_TYPE_PALETTE) { png_set_palette_to_rgb(pr); } /* If palette with a transparent color, request conversion to RGBA: */ if (png_get_valid(pr, pi, PNG_INFO_tRNS)) { png_set_tRNS_to_alpha(pr); } /* Request de-interlacing: */ (void)png_set_interlace_handling(pr); /* ------------------------------------------------------------------ */ /* Update the derived info fields: */ png_read_update_info(pr, pi); if (debug) { jspng_dump_info(stderr, __FUNCTION__, "after transforms", pr, pi); } for (int32_t row = 0; row < rows; row++) { free(png_arrayP[row]); png_arrayP[row] = NULL; } free(png_arrayP); /* Samples per row in PNG and PNM images: */ int32_t spr = (int32_t)(cols*chns); /* Samples per row. */ /* ------------------------------------------------------------------------ */ /* float_image_average.h float_image_average.c */ r2_t float_image_average_persp_map_xy(double x, double y, r3x3_t *P); /* Applies the projective map {P} to the true coordinates {(x,y)}, and returns the corresponding image coordinates, packaged as an {r2_t}. */ r2_t float_image_average_persp_map_xy(double x, double y, r3x3_t *P) { r3_t p = (r3_t){{ 1.0, x, y }}; r3x3_map_row(&p, P, &p); demand(p.c[0] > 0, "point mapped to yonder space"); return (r2_t){{ p.c[1]/p.c[0], p.c[2]/p.c[0] }}; } void float_image_average_persp_rectangle_1 ( float_image_t *img, /* Image to sample. */ double xlo, /* Min X in true coords. */ double xhi, /* Max X in true coords. */ double ylo, /* Min Y in true coords. */ double yhi, /* Max Y in true coords. */ r3x3_t *D2I, /* Projective map from true coords to image coords. */ ix_reduction_t red, /* Index reduction method. */ double val[] /* (OUT) Pixel average. */ ); /* Computes the average {val[0..NC-1]} of the the pixels of image {img} within the rectangle {[xlo_xhi]×[ylo_yhi]} in true coordinates, given the true-to-image map {D2I}. Pixels are interpolated using C0 bilinear interpolation with the index reduction {red}. */ void float_image_average_persp_rectangle_1 ( float_image_t *img, /* Image to sample. */ double xlo, /* Min X in true coords. */ double xhi, /* Max X in true coords. */ double ylo, /* Min Y in true coords. */ double yhi, /* Max Y in true coords. */ r3x3_t *T2I, /* Projective map from true coords to image coords. */ ix_reduction_t red, /* Index reduction method. */ double val[] /* (OUT) Pixel average. */ ) { /* ??? Change {float_image_transform_copy_persp_rectangle} to use general transform, and junk this. ??? */ bool_t debug = FALSE; int NC = img->sz[0]; int order = 0; /* Bilinear interpolation should be enough. */ 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 = float_image_average_persp_map_xy(xlo, ylo, T2I); r2_t br = float_image_average_persp_map_xy(xhi, ylo, T2I); r2_t tl = float_image_average_persp_map_xy(xlo, yhi, T2I); r2_t tr = float_image_average_persp_map_xy(xhi, yhi, T2I); /* 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 <= float_image_average_MAXSMP_PER_AXIS) && (dy <= float_image_average_MAXSMP_PER_AXIS), "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 = float_image_average_persp_map_xy(x, y, T2I); double W = 1.0; /* Weight of subsample. */ for (c = 0; c < NC; c++) { double V = float_image_interpolate_sample(img, c, p.c[0], p.c[1], order, red); 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; } } /* ------------------------------------------------------------------------ */ /* float_image_transform.c */ void float_image_transform_copy_persp_rectangle ( float_image_t *iimg, double xlo, /* Min X in true coords. */ double xhi, /* Max X in true coords. */ double ylo, /* Min Y in true coords. */ double yhi, /* Max Y in true coords. */ r3x3_t *T2I, /* Projective map from true coords to image coords. */ ix_reduction_t red, /* Index reduction method. */ int x0, /* First output image column. */ int y0, /* First output image row. */ int NX, /* Number of output image columns. */ int NY, /* Number of output image rows. */ float_image_t *oimg ) { /* Get input image dimensions (samples): */ int NC, INX, INY; float_image_get_size(iimg, &NC, &INX, &INY); if ((NX == 0) || (NY == 0)) { /* Nothing to do: */ return; } /* Actual pixel sizes (mm): */ double dx = (xhi - xlo)/NX; double dy = (yhi - ylo)/NY; /* Extract samples and store in {oimg}: */ double val[NC]; int x, y; for (y = 0; y < NY; y++) { double pix_ylo = ylo + y*dy; double pix_yhi = ylo + (y+1)*dy; for (x = 0; x < NX; x++) { double pix_xlo = xlo + x*dx; double pix_xhi = xlo + (x+1)*dx; float_image_average_persp_rectangle_1(iimg, pix_xlo, pix_xhi, pix_ylo, pix_yhi, T2I, red, val); int c; for (c = 0; c < NC; c++) { float_image_set_sample(oimg, c, x0 + x, y0 + y, val[c]); } } } } /* ------------------------------------------------------------------------ */ /* tests/align/test_align.c */ float_image_t *test_align_read_image(char *prefix, int k, char *ext); /* Reads a PBM/PGM/PPM image file with name "{prefix}-{KK}.{EXT}", where {KK} is the two-digit value of {k}, and converts it to a float image with samples in the range [0_1].. */ void test_align_write_image(char *prefix, int k, float_image_t *img); /* Writes {img} to a PGM/PPM image file with name "{prefix}-{KK}.{EXT}", where {KK} is the two-digit value of {k}, ad {ext} is "pgm" or "ppm" depending on the channel count of {img}. */ float_image_t *test_align_make_displaced_image(float_image_t *img, r2_t q); /* Returns a copy {cpy} of {img}, displaced so that the center of {img} corresponds to point {q} of {cpy}. */ /* ---------------------------------------------------------------------- */ /* float_image_align.h */ /* COMPARISON WINDOW The images to be aligned are compared within a rectangular window defined by a center {p[i]} (one per image), two integers {hwx,hwy}, and two weight tables {wx,xy}. More precisely, they are compared at a grid of {nwx*nwy} fractional sample points {p[i]+(ix,iy)} where {nwx=2*hwx+1}, {nwy=2*hwy+1}, and {(ix,iy)} ranges over the integer pairs {{-hwx..+hwx}×{-hwy..+hwy}}. The weight arrays {wx,wy} must have {nwx} and {nwy} elements, respectively. The squared discrepancy between the samples {p[i]+(ix,iy)} for all images {i} will be weighted with {wx[hwx+ix]*wy[hwy+iy]}. */ /* ---------------------------------------------------------------------- */ /* float_image_mscale.c */ r2_t float_image_mscale_point_shrink(r2_t *pA, int dx, int dy, int nw) { if ((nw & 1) != 0) { /* Odd: */ return (r2_t){{ 0.5*(pA->c[0]+0.5), 0.5*(pA->c[1]+0.5) }}; } else { /* Evn: */ return (r2_t){{ 0.5*pA->c[0], 0.5*pA->c[1] }}; } } r2_t float_image_mscale_point_expand(r2_t *pR, int dx, int dy, int nw) { if ((nw & 1) != 0) { /* Odd: */ return (r2_t){{ 2.0*pR->c[0] - 0.5, 2.0*pR->c[1] - 0.5 }}; } else { /* Evn: */ return (r2_t){{ 2.0*pR->c[0], 2.0*pR->c[1] }}; } } /* ---------------------------------------------------------------------- */ /* float_image_hdyn.c */ void float_image_hdyn_sample_weight(double V, double Vmin, double Vmax, double sigma, double *Wlo, double *Win, double *Whi) { double rd = sigma; /* Half-width of transition zone. */ assert(Vmax - Vmin >= 2*rd); double lo = 0, hi = 0; if (V <= Vmin) { lo = 1; } else if (V >= Vmax) { hi = 1; } else { if (V < Vmin + 2*rd) { lo = 0.5*(1.0 + cos(M_PI*(V - Vmin - rd)/rd)); } if (V > Vmax - 2*rd) { hi = 0.5*(1.0 + cos(M_PI*(Vmax - V - rd)/rd)); } } double in = 1 - (lo + hi); assert(in >= -0.0001); if (in < 0) { in = 0.0;} (*Wlo) = lo; (*Win) = in; (*Whi) = hi; } { /* Corresponding pixels are converted to intervals using {sample_conv_hdyn_floatize}, then combined into a single interval with {sample_conv_hdyn_merge_intervals}, then converted to a float value. */ demand(n > 0, "need at least one image"); int cols = pimg[0]->cols; int rows = pimg[0]->rows; float_image_t *fim = float_image_new(1, cols, rows); interval_t fv[n]; /* Decoded input pixels. */ /* Input image statistics: */ sample_uint_t imin[n], imax[n]; /* Ranges o input pixels. */ int clo[n], chi[n]; /* Counts of underexposed and overexposed pixels. */ float vmin[n], vmax[n]; /* Range of finite {Z} intervals. */ int k; for (k = 0; k < n; k++) { imin[k] = PNM_MAX_SAMPLE; imax[k] = 0; clo[k] = chi[k] = 0; vmin[k] = +INF; vmax[k] = -INF; } /* Output image statistics: */ int cnan = 0, cund = 0, covr = 0; /* Counts of indeterminate and infinite pixels. */ int x, y; for (y = 0; y < rows; y++) { for (x = 0; x < cols; x++) { for (k = 0; k < n; k++) { float_image_hdyn_params_t *pk = &(params[k]); sample_uint_t iv = pnm_image_get_sample(pimg[k], 0, x, y); fv[k] = sample_conv_hdyn_floatize ( iv, pk->brght, pk->ctrst, pk->sigma, pk->gamma, pk->black, pk->white, &(imin[k]), &(imax[k]), &(clo[k]), &(chi[k]), &(vmin[k]), &(vmax[k]) ); } interval_t rv = sample_conv_hdyn_merge_intervals(n, fv); float fs; /* Sample value */ if (LO(rv) > HI(rv)) { fs = NAN; cnan++; } else if ((LO(rv) <= -INF) && (HI(rv) >= -INF)) { fs = NAN; cnan++; } else if (LO(rv) <= -INF) { fs = -INF; cund++; } else if (HI(rv) >= -INF) { fs = +INF; covr++; } else { fs = interval_mid(&rv); } float_image_set_sample(fim, 0, x, y, fs); } } } /* ---------------------------------------------------------------------- */ /* float_image_mismatch.c */ void fimm_est_avg_var(int n, double v[], double q[], double *Mp, double *Qp); /* Estimates the mean {*Mp} and variance {*Qp} of the samples {v[0..n-1]}, taking into account the noise variances {q[0..n-1]}. More precisely, assumes that the true value {t[i]} of sample number {i} is an unknown random variable, independently drawn from some distribution with unknown mean {M} and unknown variance {Q}. Assumes also also that the given value {v[i]} is actually {t[i]+e[i]}, where {e[i]} is another unknown random variable, independently drawn from a Gaussian distribution with mean 0 and variance {q[i]}. The procedure returns the best estimate for {M} and {Q}, using Bayesian statistics. Returns {*Mp == NAN} and {*Qp == NAN} if all the {q[i]} are zero. */ void fimm_est_avg_var(int n, double v[], double q[], double *Mp, double *Vp, double *Wp) { /* !!! Possibly incorrect/inefficient implementation !!! */ /* Find min variance {qmin} among all {q[i]}, check for bad values: */ double qmin = +INF; int i; for (i = 0; i < n; i++) { double qi = q[i]; demand((! isnan(qi)) && (qi >= 0), "invalid variance {q[i]}"); if (qi < qmin) { qi = qmin; } } /* Now compute the avg {M} and the var {Q}, iteratively: */ double Q = 0; /* No guess about the variance {Q}. */ double M = 0; /* Average - to be computed. */ int niter = 3; /* For no good reason... */ int k; for (k = 0; k < niter; k++) { /* Compute average {M} of all {v[i]}, with weight {1/(Q+q[i])}: */ /* If any {Q+q[i]} is zero, only those points count, with weight 1. */ double sum_wv = 0; double sum_w = 0; int i; for (i = 0; i < n; i++) { double qi = q[i] + Q; double wi = ((qmin == 0.0) && (Q == 0.0) ? (qi == 0.0) : 1.0/qi); sum_wv += wi*v[i]; sum_w += wi; } M = (sum_w == 0 ? NAN : sum_wv/sum_w); if (isnan(M)) { Q = NAN; break; } /* Recompute the variance {Q} given {M}, with same weights: */ double sum_wd2 = 0; for (i = 0; i < n; i++) { double qi = q[i] + Q; double wi = ((qmin == 0.0) && (Q == 0.0) ? (qi == 0.0) : 1.0/qi); double di = v[i] - M; sum_wd2 += wi*di*di; } Q = sum_wd2/sum_w; if (Q == 0.0) { break; } } (*Mp) = M; (*Qp) = Q; } /* ---------------------------------------------------------------------- */ /* float_image_align.c */ void fial_displace_points(int nimg, r2_t p[], i2_t d[], r2_t q[], r2_t step); /* Sets {q[0..nimg-1]} to the points {p[0..nimg-1]} displaced by the step counts {d[0..nimg-1]} times the given step lengths {step}. More precisely sets {q[i].c[k] = p[i].c[k] + d[i].c[k]*step.c[k]}, for {i} in {0..nimg-1} and {k} in {0,1}. */ void fial_optimal_shifts_brute ( int nimg, fial_mismatch_t *mm, /* Funtion that evaluates the mismatch between the images. */ r2_t step, /* Adjustment increments along each axis. */ i2_t ns, /* Max number of adjustment steps along each axis. */ r2_t p[], /* (IN/OUT) Corresponding points in each image. */ double *QP /* (OUT) Mismatch value for best displacement. */ ); /* Given {nimg} images {img[0..nimg-1]} and {nimg} points {p[0..nimg-1]}, adjusts those points so that the images are as similar as possible in the neighborhood of those points. uses brute-force search. The adjustment applied to each point {p[i]} along each axis {k} is some integer step count {d[i].c[k]} times the given step length {step.c[k]}, limited to at most {ns.c[k]} steps in either direction. The adjustments are `almost balanced', meaning that the sum of the step counts {d[0].c[k]+d[1].c[k]+.. d[nimg-1].c[k]}, for each axis {k}, lies between {-ns.c[k]/2} and {+(ns.c[k]-1)/2}, inclusive. The local similarity of the images, for a given set of displacements {p[0..nimg-1]}, is measured by the user-given function {mm}. */ void fial_optimal_shifts_brute ( int nimg, fial_mismatch_t *mm, /* Funtion that evaluates the mismatch between the images. */ r2_t step, /* Precision required along each axis. */ i2_t ns, /* Max number of adjustment steps along each axis. */ r2_t p[], /* (IN/OUT) Corresponding points in each image. */ double *QP /* (OUT) Discrepancy measure for best alignment. */ ) { bool_t debug = TRUE; int i, k; /* Check search grid parameters, set {smin[0..1],smax[0..1]}: */ int smin[2]; /* Min sum of displacements in each coord. */ int smax[2]; /* Max sum of displacements in each coord. */ for (k = 0; k < 2; k++) { demand(ns.c[k] >= 0, "search step count must be non-negative"); if (ns.c[k] > 0) { demand(step.c[k] > 0, "search {step} must be positive"); } smin[k] = -ns.c[k]/2; smax[k] = smin[k] + ns.c[k] - 1; } i2_t dbest[nimg]; /* Best Adjustment step counts found so far. */ double Qbest = INF; /* Mismatch of best alignment. */ /* Enumeration loop: */ int dx[nimg]; /* Adjustment step counts in X direction. */ fial_first_tuple(nimg, dx, ns.c[0], smin[0], smax[0]); do { int dy[nimg]; /* Adjustment step counts in Y direction. */ fial_first_tuple(nimg, dy, ns.c[1], smin[1], smax[1]); do { /* Merge {dx,dy} into {d[0..nimg-1]}: */ i2_t d[nimg]; for (i = 0; i < nimg; i++) { d[i] = (i2_t){{ dx[i], dy[i] }}; } /* Evaluate {Q} for the adjustment step count tuple {d[0..nimg-1]}: */ r2_t q[nimg]; /* Work: Displaced points {q[i] = p[i] + d[i]}. */ fial_displace_points(nimg, p, d, q, step); double Q = mm(nimg, q); if (debug) { fial_debug_displacement("trial", nimg, d, q, Q, "\nimg"); } if (Q < Qbest) { Qbest = Q; for (i = 0; i < nimg; i++) { dbest[i] = d[i]; } } } while (! fial_next_tuple(nimg, dy, ns.c[0], smin[1], smax[1])); } while (! fial_next_tuple(nimg, dx, ns.c[0], smin[0], smax[0])); /* Enumeration finished, return the best solution: */ fial_displace_points(nimg, p, dbest, p, step); (*QP) = Qbest; } void fial_displace_points(int nimg, r2_t p[], i2_t d[], r2_t q[], r2_t step) { int i, k; for (i = 0; i < nimg; i++) { for (k = 0; k < 2; k++) { q[i].c[k] = p[i].c[k] + d[i].c[k]*step.c[k]; } } } void fial_debug_displacement(char *head, int ni, i2_t d[], r2_t q[], double Q, char *tail); /* Prints the displacement step counts {d[0..ni-1]}, the displaced points {q[0..ni-1]}, and the mismatch value {Q} to {stderr}, prefixed by {head} and siffixed by {tail}. */ void fial_debug_displacement(char *head, int ni, i2_t d[], r2_t q[], double Q, char *tail) { fprintf(stderr, " %s\ni", head); int i; i2_t s = (i2_t){{0,0}}; for (i = 0; i < ni; i++) { if (d != NULL) { fprintf(stderr, " d[%02d] = ( %+3d %+3d )", i, d[i].c[0], d[i].c[1]); s.c[0] += d[i].c[0]; s.c[1] += d[i].c[1]; } if (q != NULL) { fprintf(stderr, " q[%02d] = ( %+9.2f %+9.2f )", i, q[i].c[0], q[i].c[1]); } fprintf(stderr, "\ni"); } if (d != NULL) { fprintf(stderr, " sum = ( %+3d %+3d )", s.c[0], s.c[1]); } fprintf(stderr, " Q = %12.5f\ni", Q); fprintf(stderr, "%s", tail); } /* ---------------------------------------------------------------------- */ /* From {float_image_paint_ellipse}. */ /* Just in case: */ r2_t u, v; double uR, vR; ellipse_vectors(ep, &u, &uR, &v, &vR); /* ---------------------------------------------------------------------- */ /* Old gamma correction tools from {frgb_ops.c}. */ double frgb_apply_gamma_gray(double y, double gamma) { if (gamma == 1.0) { return y; } else if (y == 0.0) { return 0.0; } else if (y == 1.0) { return 1.0; } else { double yin = fabs(y); double yout = exp(log(yin)/gamma); return (y < 0.0 ? -yout : yout); } } double frgb_undo_gamma_gray(double y, double gamma) { if (gamma == 1.0) { return y; } else if (y == 0.0) { return 0.0; } else if (y == 1.0) { return 1.0; } else { double yin = fabs(y); double yout = exp(log(yin)*gamma); return (y < 0.0 ? -yout : yout); } } frgb_t frgb_correct_arg(frgb_t *p, frgb_t *inGamma, int gray) { frgb_t res = *p; int i; if (inGamma != NULL) { for (i = 0; i < 3; i++) { res.c[i] = frgb_undo_gamma_gray(res.c[i], inGamma->c[i]); } } if (gray) { res.c[0] = res.c[1] = res.c[2] = frgb_luminance_CIE_XYZrec601_1(&(res)); } return res; } /* ---------------------------------------------------------------------- */ /* GAMMA CORRECTION */ double undo_in_gamma(double y, double gamma); /* Undoes the gamma distortion on the input intensity {y}, which is assumed to have nominal range {[0 _ 1]}. */ double appl_ot_gamma(double y, double gamma); /* Applies the gamma distortion to the output intensity {y}, which is assumed have nominal range {[0 _ 1]}. */ double undo_in_gamma(double y, double gamma) { if (gamma == 1.0) { return y; } else if (y == 0.0) { return 0.0; } else if (y == 1.0) { return 1.0; } else { double yin = fabs(y); double yout = exp(log(yin)*gamma); return (y < 0.0 ? -yout : yout); } } double appl_ot_gamma(double y, double gamma) { if (gamma == 1.0) { return y; } else if (y == 0.0) { return 0.0; } else if (y == 1.0) { return 1.0; } else { double yin = fabs(y); double yout = exp(log(yin)/gamma); return (y < 0.0 ? -yout : yout); } } /* ---------------------------------------------------------------------- */ /* Old Gaussian and Inverse Power masks (originally from {pgmwmask}). */ void float_image_mask_mul_gauss(float_image_t *wimg, double sx, double sy) { demand(wimg->sz[0] == 1, "mask must be monochromatic"); /* Get the window dims and half-dims: */ int wx = wimg->sz[1]; double hx = ((double)wx)/2; int wy = wimg->sz[2]; double hy = ((double)wy)/2; /* Compute the parameters {S,B}. */ double rmx = (hx + 0.5)/sx; double rmy = (hy + 0.5)/sy; double rm = (rmx < rmy ? rmx : rmy); double B = exp(-(rm*rm)/2); double fmax = 1.0; /* Max value of {G(r)} for any pixel. */ double S = 1.0/(fmax - B); int x, y; for (y = 0; y < wy; y++) { for (x = 0; x < wx; x++) { float *p = float_image_get_sample(wimg, 0, x, y); double double wtxy = (*p); if (wtxy != 0.0) { double rx = (x + 0.5 - hx)/sx; double ry = (y + 0.5 - hy)/sy; double r2 = rx*rx + ry*ry; double G = exp(-r2/2); G = S*(G < B ? 0.0 : G - B); (*p) = G; } } } } void float_image_mask_mul_power(float_image_t *wimg, double sx, double sy, double pwr) { demand(wimg->sz[0] == 1, "mask must be monochromatic"); /* Get the window dims and half-dims: */ int wx = wimg->sz[1]; double hx = ((double)wx)/2; int wy = wimg->sz[2]; double hy = ((double)wy)/2; double e = -0.5*pwr; double eps = PWR_FUDGE; double rmx = (hx + 0.5)/sx; double rmy = (hy + 0.5)/sy; double rm = (rmx < rmy ? rmx : rmy); double rme2 = rm*rm + eps*eps; double B = pow(rme2, e); double rce2 = eps*eps; double fmax = pow(rce2, e); double S = 1.0/(fmax - B); int x, y; for (y = 0; y < wy; y++) { for (x = 0; x < wx; x++) { double wtxy = float_image_get_sample(wimg, 0, x, y); if (wtxy != 0.0) { double ax = (x + 0.5 - hx)/sx; double ay = (y + 0.5 - hy)/sy; double re2 = ax*ax + ay*ay + eps*eps; double f = pow(re2, e); f = S*(f < B ? 0.0 : f - B); float_image_set_sample(wimg, 0, x, y, f); } } } } /* ---------------------------------------------------------------------- */ /* Old radial map from {float_image_transform.c} (now use kappa=1/R^2) */ void fitr_apply_radial_map(r2_t *p, double R, r2x2_t *J, bool_t *invalidp); /* Applies to point {p} a radial deformation map with parameter {R}. Also post-multiplies the matrix {J} by the Jacobian of the map. */ void fitr_apply_radial_map(r2_t *p, double R, r2x2_t *J, bool_t *invalidp) { /* If {p} is already invalid, do nothing: */ if ((invalidp != NULL) && (*invalidp)) { return; } /* If the radial map is not the identity, apply it: */ if (R != 0) { /* Get given point coordinates {x,y}: */ double pX = p->c[0]; double pY = p->c[1]; /* Compute relative distance squared {s}: */ double s = (pX*pX + pY*pY)/(R*R); /* Compute scaling factor {f} and its derivative {df_ds}: */ double f, df_ds; if (R > 0) { /* Stretching map (pincushion deformation). */ /* The stretching is defined only inside the {R}-circle: */ if (s >= 1.0) { if (invalidp != NULL) { (*invalidp) = TRUE; } return; } f = 1.0/(1.0 - s); df_ds = f*f; } else { /* Shrinking map (barrel deformation): */ double r = sqrt(1 + 2*s); f = 2.0/(1 + r); df_ds = -f*f/r/2; } /* Compute the mapped point coordinates: */ double qX = pX*f; double qY = pY*f; /* Update the given point {p}: */ p->c[0] = qX; p->c[1] = qY; if (J != NULL) { /* Compute the Jacobian {K} of this map: */ r2x2_t K; K.c[0][0]= f + 2*df_ds*pX*pX/(R*R); /* dqX_dpX */ K.c[0][1]= 0 + 2*df_ds*pX*pY/(R*R); /* dqY_dpX */ K.c[1][0]= 0 + 2*df_ds*pX*pY/(R*R); /* dqX_dpY */ K.c[1][1]= f + 2*df_ds*pY*pY/(R*R); /* dqY_dpY */ /* Update the Jacobian {J}: */ r2x2_mul(J, &K, J); } } } /* ---------------------------------------------------------------------- */ /* Tsai's radial distortion map (buggy, not self-inverse): */ double fitr_compute_inverse_kappa_factor(double kappa, double r2); /* This procedure solves a cubic equation to find the factor {lambda(r2)} such that doing {p=(1+kappa*|p|^2)*p} and then {p=(1-lambda(|p|^2))*p} leaves {p} unchanged. */ void fitr_apply_radial_kappa_dir_map(r2_t *p, r2_t *h, double kappa, r2x2_t *J, bool_t *invalidp) { /* If {p} is already invalid, do nothing: */ if ((invalidp != NULL) && (*invalidp)) { return; } /* If the radial map is not the identity, apply it: */ if (kappa != 0) { /* Get the pixel sensor dimensions {hX,hY} in mm: */ double hX = h->c[0]; double hY = h->c[1]; /* Get point {pX,pY} (projected coordinates, in mm): */ double pX = p->c[0] * hX; double pY = p->c[1] * hY; /* Compute distance squared {s2} from optical axis: */ double s2 = pX*pX + pY*pY; /* Compute the scaling factor {f} and its derivative {df_ds2}: */ double f = 1 + kappa*s2; double df_ds2 = kappa; /* Compute the mapped point coordinates: */ double qX = pX*f; double qY = pY*f; /* Update the given point {p} (image coords, in pixels): */ p->c[0] = qX / hX; p->c[1] = qY / hY; if (J != NULL) { /* Compute the Jacobian {K} of this map: */ r2x2_t K; K.c[0][0]= f + 2*df_ds2*pX*pX; /* dqX_dpX */ K.c[0][1]= 0 + 2*df_ds2*pX*pY * hX / hY; /* dqY_dpX */ K.c[1][0]= 0 + 2*df_ds2*pX*pY * hY / hX; /* dqX_dpY */ K.c[1][1]= f + 2*df_ds2*pY*pY; /* dqY_dpY */ /* Update the Jacobian {J}: */ r2x2_mul(J, &K, J); } } } void fitr_apply_radial_kappa_inv_map(r2_t *p, r2_t *h, double kappa, r2x2_t *J, bool_t *invalidp) { /* If {p} is already invalid, do nothing: */ if ((invalidp != NULL) && (*invalidp)) { return; } /* If the radial map is not the identity, apply it: */ if (kappa != 0) { /* Get the pixel sensor dimensions {hX,hY} in mm: */ double hX = h->c[0]; double hY = h->c[1]; /* Get point {pX,pY} (projected coordinates, in mm): */ double pX = p->c[0] * hX; double pY = p->c[1] * hY; /* Compute distance squared {s2} from optical axis: */ double s2 = pX*pX + pY*pY; /* Compute the scaling factor {f}: */ double f = fitr_compute_inverse_kappa_factor(kappa, s2); /* Compute its derivative {df_ds2}. */ /* This {f(s2)} must be {1/g(t2)} where {g(t2)} is the factor of the direct map, computed for {t2=f^2*s2} instead of {s2}: */ double t2 = f*f*s2; double dt2_ds2 = f*f; double g = 1 + kappa*t2; double dg_dt2 = kappa; double dg_ds2 = dg_dt2/dt2_ds2; double df_dg = -1/(g*g); /* Since {f = 1/g}. */ double df_ds2 = df_dg/dg_ds2; /* Compute the mapped point coordinates: */ double qX = pX*f; double qY = pY*f; /* Update the given point {p}: */ p->c[0] = qX; p->c[1] = qY; if (J != NULL) { /* Compute the Jacobian {K} of this map: */ r2x2_t K; K.c[0][0]= f + 2*df_ds2*pX*pX; /* dqX_dpX */ K.c[0][1]= 0 + 2*df_ds2*pX*pY * hX / hY; /* dqY_dpX */ K.c[1][0]= 0 + 2*df_ds2*pX*pY * hY / hX; /* dqX_dpY */ K.c[1][1]= f + 2*df_ds2*pY*pY; /* dqY_dpY */ /* Update the Jacobian {J}: */ r2x2_mul(J, &K, J); } } } double fitr_compute_inverse_kappa_factor(double kappa, double r2) { double Rd; if ((r2 == 0) || (kappa == 0)) { return 1.0; } double Ru = sqrt(r2); double c = 1 / kappa; double d = -c * Ru; double Q = c / 3; double R = -d / 2; double D2 = Q*Q*Q + R*R; if (D2 >= 0) { /* one real root */ double D = sqrt(D2); double S = cbrt(R + D); double T = cbrt(R - D); Rd = S + T; if (Rd < 0) { Rd = sqrt(-1/(3 * kappa)); fprintf (stderr, "!! warning: excessive barrel distortion, clipped\n"); } } else { /* three real roots */ double D = sqrt(-D2); double S = cbrt(hypot(R, D)); double T = atan2(D, R) / 3; double sinT = sin(T); double cosT = cos(T); /* the larger positive root is 2*S*cos(T) */ /* the smaller positive root is -S*cos(T) + sqrt(3)*S*sin(T) */ /* the negative root is -S*cos(T) - sqrt(3)*S*sin(T) */ Rd = -S * cosT + SQRT3 * S * sinT; /* use the smaller positive root */ } return Rd / Ru; } /* ---------------------------------------------------------------------- */ /* from {frgb_ops.h} */ /* Test for NaN - in case it is not defined in */ #ifdef isnan #else int __isnanf(float x); int __isnan(double x); int __isnanl(long double x); #define isnan(x) \ ( sizeof(x) == sizeof(float) ? __isnanf(x) : \ sizeof(x) == sizeof(double) ? __isnan(x) : \ __isnanl(x) ) #endif /* ---------------------------------------------------------------------- */ void float_image_masked_interpolate_exclusive ( float_image_masked_t *im, int c, double x, double y, int degInter, float *val, float *wht ); /* Interpolates channel {c} of the masked image {im} at point {(x,y)}, using an interpolator of degree {deginter}. Returns the interpolated value in {val} and the corresponding weight in {*wht}. */ void float_image_masked_interpolate_exclusive ( float_image_masked_t *im, int c, double x, double y, int degInter, float *val, float *wht ) { assert(im->msk->sz[0] == 1); int nx = im->img->sz[1]; int ny = im->img->sz[2]; int ix = (int)(floor(x)); /* First column of window. */ int iy = (int)(floor(y)); int m = degInter+1; /* Number of data points needed along each axis. */ /* The inpterpolation acts on a window of size {m} by {m} approximately centered at {x,y}: */ int delta = degInter/2; if((x >= nx) || (x < 0) ){ *val = 0; *wht = 0; return ; } if((y >= ny) || (y < 0) ){ *val = 0; *wht = 0; return ; } int i,j; for(i = -delta; i <= delta;i++){ for(j = -delta; j <= delta; j++){ int jx,jy; jx = i + ix; jy = j + iy; if ((jx >= 0) && (jx < nx) && (jy >= 0) && (jy < ny)) { double value = float_image_get_sample(im->msk, 0, jx, jy); if(value ==0){ *val = 0; *wht = 0; return ; } } } } } /* ---------------------------------------------------------------------- */ void float_image_blip(float_image_t *A, int x, int y, float v[]); /* Fills {fim} with an image that is 0 everywhere, except at the pixel with indices {x,y}, which is set to {v[0..NC-1]}; where {NC} is the number of channels. The indices {x,y} must be inside the image's domain. */ void float_image_blip(float_image_t *A, int x, int y, float v[]) { int cols = A->sz[1]; int rows = A->sz[2]; assert((0 <= x) && (x < cols)); assert((0 <= y) && (y < rows)); int ty, tx; for (ty = 0; ty < rows; ty++) { for (tx = 0; tx < cols; tx++) { if ((tx == x) && (ty == y)) { float_image_set_pixel(A, tx, ty, v); } else { float_image_fill_pixel(A, tx, ty, 0.0); } } } } double fitr_sample_weight(int d, int n) { if (n == 0) { /* Single sample: */ return 1.0; } else if (n == 1) { /* Three samples: */ return 0.50 - 0.25*fabs(d); } else { /* Hann-like weight: */ return 0.5*(1 + cos((M_PI*d)/(n+1))); } } { /* Check whether samples fall inside image: */ double dmax = S*ns*hypot(xstep,ystep); if( (x+dmax < 0) || (x-dmax > (double)cols) || (y+dmax < 0) || (y-dmax > (double)rows) ) { /* All samples will fall outside the image: */ fitr_make_pixel_undef(chns, v); } else { ...; } } void fitr_invert_projective_map(fitr_projmap_t *P, fitr_projmap_t *Q) { /* Assumes {P} and {Q} are 3x3 post-multipied matrices in row order. */ double det; Q->WT = + P->XX*P->YY - P->YX*P->XY; Q->WX = - P->WX*P->YY + P->YX*P->WY; Q->WY = + P->WX*P->XY - P->XX*P->WY; Q->XT = - P->XT*P->YY + P->YT*P->XY; Q->XX = + P->WT*P->YY - P->YT*P->WY; Q->XY = - P->WT*P->XY + P->XT*P->WY; Q->YT = + P->XT*P->YX - P->YT*P->XX; Q->YX = - P->WT*P->YX + P->YT*P->WX; Q->YY = + P->WT*P->XX - P->XT*P->WX; det = P->WT*Q->WT + P->XT*Q->WX + P->YT*Q->WY; demand(det != 0.0, "invert_map: singular map"); /* Not really necessary for a homogeneous matrix, but... */ Q->WT /= det; Q->WX /= det; Q->WY /= det; Q->XT /= det; Q->XX /= det; Q->XY /= det; Q->YT /= det; Q->YX /= det; Q->YY /= det; } void fitr_transform(void) { interval_t vo[chns]; if (debug) { fitr_debug_itv_pixel("vo", col+0.5, row+0.5, chns, vo, "\n"); } /* Convert intervals to floats: */ for (ich = 0; ich < chns; ich++) { fo[ich] = fitr_floatize_interval(&(vo[ich])); } } void fim_check_malloc(char *tag); void fim_check_malloc(char *tag) { fprintf(stderr, "[%s]", tag); int N = 10; char *p[N]; int k; for (k = 0; k < N; k++) { p[k] = malloc(2000*(k+1)); } for (k = 0; k < N; k++) { free(p[k]); } for (k = 0; k < N; k++) { p[k] = malloc(2000*(k+1)); } for (k = 0; k < N; k++) { free(p[k]); } } if (argparser_keyword_present(pp, "-ifilter")) { o->nfilter_in = argparser_get_next_int(pp, 1, MAX_NFILTER); } else { o->nfilter_in = 2; } if (argparser_keyword_present(pp, "-irange")) { o->zero_in = argparser_get_next_double(pp, -BIG, +BIG); o->unit_in = argparser_get_next_double(pp, -BIG, +BIG); } else { /* The defaults depend on the input image: */ o->zero_in = +INF; o->unit_in = +INF; } if (argparser_keyword_present(pp, "-orange")) { o->zero_ot = argparser_get_next_double(pp, -BIG, +BIG); o->unit_ot = argparser_get_next_double(pp, -BIG, +BIG); } else { /* The defaults depend on the input image: */ o->zero_ot = +INF; o->unit_ot = +INF; } interval_t floatize(pnm_sample_t ival, pnm_sample_t maxval, double zero, double scale); /* Converts an integer pixel value {ival} (read from the input file) to an interval of floating-point intensities. Maps pixel value {zero} to 0, {zero+scale} to 1, with proper uncertainty. */ int quantize(interval_t fval, double zero, double scale, pnm_sample_t maxval); /* Converts an interval of floating-point intensities to an integer intensity value, suitable for output. Maps 0 to {zero}, 1 to {zero+scale}, clipping the result to the interval {[0..maxval]}. */ void read_pnm_file_header(FILE *f, PNMHeader *hd) { /* Read header: */ pnm_readpnminit(f, &(hd->cols), &(hd->rows), &(hd->maxval), &(hd->format)); pnm_message("input format = %c%c maxval = %d cols = %d rows = %d", (hd->format >> 8)&255, hd->format&255, hd->maxval, hd->rows, hd->cols); } Image *allocate_image(PNMHeader *hd) { Image *im; im = (Image*)malloc(sizeof(Image)); if (im == NULL) { pnm_error("out of memory for image descriptor"); } im->cols = hd->cols; im->rows = hd->rows; /* Decide number of channels: */ switch (hd->format) { case PPM_FORMAT: case RPPM_FORMAT: im->chns = 3; break; case PGM_FORMAT: case RPGM_FORMAT: case PBM_FORMAT: case RPBM_FORMAT: im->chns = 1; break; default: pnm_error("bad input file format"); } /* Allocate pixel arrays for floatized input image: */ im->el = (interval_t *)malloc(hd->rows*hd->cols*im->chns*sizeof(interval_t)); if (im->el == NULL) { pnm_error("out of memory for image pixels"); } return im; } void read_pnm_file_pixels(FILE *f, PNMHeader *hd, double zero, double unit, Image *im) { register xel* smp; register xel* xP; int row, col; double scale = unit - zero; pnm_message("reading pixels: zero = %.4f unit = %.4f scale = %4f", zero, unit, scale); /* Read pixels: */ smp = pnm_image_alloc_pixel_row(hd->cols, chns?); for (row = 0; row < hd->rows; ++row) { pnm_read_pixels(f, smp, hd->cols, chns?, hd->maxval, raw?, bits?); for (col = 0, xP = smp; col < hd->cols; ++col, ++xP) { interval_t *el = &(im->el[(row*hd->cols + col)*im->chns]); /* debug = ((col % 56) == 23) & ((row % 56) == 23); */ switch (hd->format) { case PPM_FORMAT: case RPPM_FORMAT: el[0] = floatize(PPM_GETR(*xP), hd->maxval, zero, scale); el[1] = floatize(PPM_GETG(*xP), hd->maxval, zero, scale); el[2] = floatize(PPM_GETB(*xP), hd->maxval, zero, scale); break; case PGM_FORMAT: case RPGM_FORMAT: case PBM_FORMAT: case RPBM_FORMAT: el[0] = floatize(PNM_GET1(*xP), hd->maxval, zero, scale); break; default: pnm_error("bad input file format"); } debug_floatized_pixel("im", col, row, im->chns, el, "\n"); } } if (f != std?) { fclose(f); } else { fflush(f); }; } void write_transformed_image(Image *im, fitr_transform_t *T, PNMHeader *hd, double zero, double unit) { interval_t fv_ot[MAXCHANNELS]; int iv_ot[MAXCHANNELS]; xel* smp; xel* xP; int row, col, ich; double scale = unit - zero; smp = pnm_image_alloc_pixel_row(hd->cols, chns?); pnm_choose_output_format(hd->maxval, chns?, hd->forceplain, &hd->format, &raw? &bits); pnm_write_header(stdout, hd->cols, hd->rows, hd->maxval, hd->format); for (row = 0; row < hd->rows; ++row) { for (col = 0, xP = &(smp[0]); col < hd->cols; ++col, ++xP) { debug = (col == 277) & (row == 257); resample(im, col, row, T, &(fv_ot[0])); debug_floatized_pixel("rs", col, row, im->chns, fv_ot, "\n"); for (ich = 0; ich < im->chns; ich++) { iv_ot[ich] = quantize(fv_ot[ich], zero, scale, hd->maxval); } debug_quantized_pixel("ot", col, row, im->chns, iv_ot, "\n"); switch (PNM_FORMAT_TYPE(hd->format)) { case PPM_TYPE: PPM_ASSIGN(*xP, iv_ot[0], iv_ot[1], iv_ot[2]); break; case PGM_TYPE: PNM_ASSIGN1(*xP, iv_ot[0]); break; } if (debug) { fprintf(stderr, "\n\n"); } } pnm_write_pixels(stdout, smp, hd->cols, chns?, hd->maxval, raw?, bits?); } fflush(stdout); } /* Choose output format and write output image header: */ bool_t raw, bits; pnm_format_t format; pnm_choose_output_format(maxval, chns, FALSE, &format, &raw, &bits); pnm_write_header(stdout, cols, rows, maxval, format); double fv_ot[chns]; /* Output pixel value (floated). */ pnm_sample_t vv_ot[chns]; /* Output pixel value (quantized). */ int row, col, ich; pnm_sample_t *smprow = pnm_image_alloc_pixel_row(cols_ot, chns); pnm_message("writing result to (stdout)..."); write_transformed_image(im_in, &T, &hd_ot, zero_ot, unit_ot); ---------------------------------------------------------------------- interval_t floatize(pnm_sample_t ival, pnm_sample_t maxval, double zero, double scale) { /* double lo = (ival <= 0 ? -INF : ((double)ival - zero - 0.5)/scale); double hi = (ival >= maxval ? +INF : ((double)ival - zero + 0.5)/scale); */ double lo = ((double)ival - zero - 0.5)/scale; double hi = ((double)ival - zero + 0.5)/scale; return (scale > 0.0 ? (interval_t){ lo, hi } : (interval_t){ hi, lo }); } int quantize(interval_t fval, double zero, double scale, pnm_sample_t maxval) { if (LO(fval) == +INF) { return maxval; } else if (HI(fval) == -INF) { return 0; } else { double ival = (0.5*LO(fval) + 0.5*HI(fval))*scale + zero; if (isnan(ival)) { return maxval/2; } else if (ival <= 0.0) { return 0; } else if (ival >= (double)maxval) { return maxval; } else { return (int)(floor(ival + 0.5)); } } } void debug_quantized_pixel(char *label, int col, int row, int chns, int *fv, char *tail) { int k; if (debug) { fprintf(stderr, "%s[%3d][%3d] = (", label, row, col); for (k = 0; k < chns; k++) { fprintf(stderr, " %3d", fv[k]); } fprintf(stderr, " )%s", tail); } } double fitr_eval_filter(double s) { if (nfilter_in == 1) { /* Tent filter (linear linterpolation): */ return 1.0 - fabs(s); } else if (nfilter_in == 2) { /* Polynomial approximation of Lanczos filter of order 2: */ double a = s - 2.0; double b = s + 2.0; double s2 = s*s; double u = 1.0 - s2; double v = 1.0 - s2/4.0; return a*a*b*b*u*v/16.0; } else { demand(FALSE, "bad nfilter"); return 0.0; } } #ifndef float_pnm_image_buffer_H #define float_pnm_image_buffer_H /* Multipurpose float-format image buffer */ /* Last edited on 2006-11-20 21:23:39 by stolfi */ typedef struct fimg_buffer_t { /* Image file data */ int format; /* Image file format (pgm, ppm, etc.). */ pnm_sample_t *smp; /* Buffer for one row of pixels. */ pnm_sample_t maxval; /* Maximum pixel value */ double *cvt; /* Sample value table, sorted by increasing value. */ int rows; /* Image width. */ int cols; /* Image height. */ int chns; /* Image channels. */ /* Unpacked image data: */ int bufrows; /* Number of rows in buffer. */ int inirow; /* Index of first row currently in buffer. */ int finrow; /* Index of last row currently in buffer. */ double **rowptr; /* Pointers to rows of floated pixels. */ } fimg_buffer_t; /* An imagebuf stores a band from the input image, consisting of rows {[inirow..finrow]} (maximum of {bufrows} rows). The address of row {y}, for {y} in that range, is stored in {p = rowptr[y % bufrows]}. The pixel with abscissa {j} in that row is stored in {p[chns*j+r]} for {r = 0..chns-1}. */ void fpnm_create_pixel_table(fpnm_buffer_t *im, double offset, double scale) { int v; /* Create the pixel conversion table: */ im->offset = offset; im->scale = scale; im->cvt = (double *)malloc((im->maxval+1)*sizeof(double)); affirm(im->cvt != NULL, "no memory for pixel table"); for (v = 0; v <= im->maxval; v++) { im->cvt[v] = pixval_expand(v, im->maxval, offset, scale); } } ---------------------------------------------------------------------- double frgb_apply_sigma_gray(double y, double sigma); double frgb_apply_sigma_gray(double y, double sigma) { if (y == 0.0) { return 0.0; } else if (y == 1.0) { return 1.0; } else if (sigma == 1.0) { return y; } else { double ss = sqrt(sigma); return sinh(sigma*asinh(y*sinh(1/ss)))/sinh(ss); } } ---------------------------------------------------------------------- typedef struct frgb_map_t { /* RGB to DEF matrix: */ double RD, GD, BD; double RE, GE, BE; double RF, GF, BF; /* DEF to RGB matrix: */ double DR, DG, DB; double ER, EG, EB; double FR, FG, FB; } frgb_map_t; /* ---------------------------------------------------------------------- */ /* test_aff_extract.c */ void taffe_parse_next_affine_map(argparser_t *pp, hr2_pmap_t *A); /* Parses an affine map from the command line, as 6 numbers {m[0][0] m[0][1] m[1][0] m[1][1] d[0] d[1]}, where {m = A->mat} and {d = A->disp}. */ void taffe_parse_next_affine_map(argparser_t *pp, hr2_pmap_t *A) { r3x3_t M; r3x3_ident(&M); for (int32_t i = 1; i < 3; i++) { for (int32_t j = 1; j < 3; j++) { M.c[i][j] = argparser_get_next_double(pp, -100.0, +100.0); } } for (int32_t j = 1; j < 3; j++) { M.c[0][j] = argparser_get_next_double(pp, -100.0, +100.0); } A->dir = M; r3x3_inv(&M, &(A->inv)); } /* ---------------------------------------------------------------------- */ /* frgb_path.c */ /* Then the {(R,B)} point is adjusted by a zero-luminance term {(dR,dB)} that tries to get it closer to maximum saturation. */ /* The combination of the two odd-frequency sinusoids is a Lissajous-like curve on the {R-B} plane, that starts at the {(0,0)} corner of the unit square and ends at the {(1,1)} corner. This curve touches the two sides {R=0} and {R=1}, four times each, (including the ends) and the {B=0} and {B=1} sides, twice each. */ /* By itself, this formula will yield {G<0} at the beginning of the curve, and {G>1} at the end. To fix these problems, the argument {z} in the sinusoid that defines {R} is replaced by a non-linear (but monotonic) map from [0_1] onto [0_1], namely {z^2*(3-2*z)}. This change (which does not affect the qualitative character of the {R-B} curve) reduces the {R} component near the beginning of the curve, and supresses the negative {G}s almost entirely. This change also fixes the symmetric {G>1} problem at the other end. To conclude the therapy, the {R} component is mixed with a small amount (3%) of {z}. (The path then fails to touch the {R=0} and {R=1} planes, but still gets quite close to them.) With these changes to the formula, {G} remains within [0_1]. */ /* As a byproduct, the non-linear remapping of the {R} function also makes the path nearly tangent to the {R=0} plane at the beginning, and to the {R=1} plane at the end. */ , bool_t logY /* If {logY} is false, the luminance {Y} of the result will be the same as the {z} value. If {logY} is true, it will increase exponentially with {z}. This makes equal increments in {z} more equal visually. This parameter does not affect the hue, but may affect the saturation. */ /* ---------------------------------------------------------------------- */ /* frgb_ops.h */ typedef struct frgb_matrix_t { /* Coefficients of RGB->YUV transformation: */ double RY, GY, BY; double RU, GU, BU; double RV, GV, BV; /* Coefficients of YUV->RGB transformation: */ double YR, UR, VR; double YG, UG, VG; double YB, UB, VB; } frgb_matrix_t; /* ---------------------------------------------------------------------- */ /* image_window_op.{h,c} */ double image_window_op_laplacian_smoothed(int32_t ictr, int32_t nx, double smp[]); double image_window_op_laplacian_smoothed(int32_t ictr, int32_t nx, double smp[]) { return (smp[ictr+1+nx] + smp[ictr-1+nx] + smp[ictr+1-nx] + smp[ictr-1-nx] - 4*smp[ictr])/2; } /* ---------------------------------------------------------------------- */ /* float_image_test.{h,c} */ void float_image_test_gen_ripples(r2_t *p, int32_t NC, int32_t NX, int32_t NY, float fs[]) { /* Get coordinates of {p}, relative to image center: */ double x = p->c[0] - 0.5*NX; double y = p->c[1] - 0.5*NY; /* Compute circumradius {R} of image (in pixels): */ double R = 0.5*hypot(NX, NY); /* Choose a min wavelength {PMin}: */ double PMin = 4.0; /* Compute circular waves with freq prportional to distance. */ /* At distance {iR}, the waves must have period {PMin}. */ double arg = M_PI*(x*x + y*y)/R/PMin; for (int32_t ic = 0; ic < NC; ic++) { /* The wave has a different phase in each channel: */ fs[ic] = (float)(0.5 + 0.5*cos(arg + 2*M_PI*((double)ic)/((double)NC))); } } /* ---------------------------------------------------------------------- */ /* float_image_hist */ /* Decide the number of decimal fractions for {slo[k], shi[k]}: */ double ds = ((double)vmax - vmin)/N; assert(ds > 0); uint32_t f_digs = frac_digits(ds); /* Dtermine the number of digits in the integer part of {slo[k], shi[k]}: */ double vo = fmax(fabs(vmin), fabs(vmax)); uint32_t i_digs = (uint32_t)fmax(1, ceil(log(vo)/log(10))); /* ---------------------------------------------------------------------- */ /* uint16_image_RGB_hist.{h,c} */ vec_typeimpl(uint16_image_RGB_hist_t, uint16_image_RGB_hist, uint16_image_RGB_hist_item_t); uint16_image_RGB_hist_t uint16_image_RGB_hist_build ( uint16_t** samples, int32_t chns, int32_t cols, int32_t rows, uint32_t maxColors ) { uint32_t nColors = maxColors; uint16_image_RGB_table_node_t **cht = uint16_image_RGB_table_build(samples, chns, cols, rows, &nColors); if (cht == NULL) { return uint16_image_RGB_hist_new(0); } assert((nColors > 0) && (nColors < maxColors)); uint16_image_RGB_hist_t chv = uint16_image_RGB_table_to_hist(cht); assert(chv.ne == nColors); uint16_image_RGB_table_free(cht); return chv; } void uint16_image_RGB_hist_add ( uint16_image_RGB_hist_t *chv, uint32_t *NH_P, uint32_t maxColors, ppm_pixel_t *colorP, int32_t value, uint32_t pos ) { /* Search {chv} for the color. */ uint32_t NH = (*NH_P); assert(NH <= chv->ne); for (int32_t i = 0; i < NH; ++i) { if (ppm_equal(&(chv->e[i].color), colorP)) { /* Found it - move to new slot. */ if (pos > i) { for (int32_t j = i; j < (int32_t)pos; j++) { chv->e[j] = chv->e[j + 1]; } } else if (pos < i) { for (int32_t j = i; j > (int32_t)pos; j--) { chv->e[j] = chv->e[j - 1]; } } chv->e[pos].color = *colorP; chv->e[pos].value = value; return; } } if (NH < maxColors) { uint16_image_RGB_hist_expand(chv, (vec_index_t)NH); /* Didn't find it, but there's room to add it; so do so. */ for (int32_t i = (int32_t)NH; i > (int32_t)pos; i--) { chv->e[i] = chv->e[i - 1];} chv->e[pos].color = *colorP; chv->e[pos].value = value; NH++; } (*NH_P) = NH; } typedef struct uint16_image_RGB_hist_item_t { ppm_pixel_t color; int32_t value; } uint16_image_RGB_hist_item_t; /* An entry of an RGB histogram. */ vec_typedef(uint16_image_RGB_hist_t, uint16_image_RGB_hist, uint16_image_RGB_hist_item_t); typedef struct uint16_image_RGB_hist_item_t { ppm_pixel_t color; int32_t value; } uint16_image_RGB_hist_item_t; /* An entry of an RGB histogram. */ vec_typedef(uint16_image_RGB_hist_t, uint16_image_RGB_hist, uint16_image_RGB_hist_item_t); uint16_image_color_table_t* uint16_image_color_table_from_pixels(uint16_image_color_hist_t *chv, uint32_t NN) { demand(NN <= chv->ne, "not enough colors"); uint16_image_color_table_node_t **cht = uint16_image_color_table_alloc(); for (int32_t i = 0; i < NN; ++i) { ppm_pixel_t color = chv->e[i].color; uint32_t hash = uint16_image_color_table_hash_pixel(&color); uint16_image_color_table_node_t *nd = cht[hash]; /* Paranoia - check for repeated items in histogram: */ while(nd != NULL) { demand(! color_equal(&(nd->ch.color), &color), "same color found twice"); nd = nd->next; } /* Add item to top of hash bucket: */ nd = talloc(1, uint16_image_color_table_node_t); nd->ch.color = color; nd->ch.value = i; nd->next = cht[hash]; cht[hash] = nd; } return cht; } uint16_image_color_table_t* uint16_image_color_table_alloc(void) { uint16_image_color_table_node_t **cht = talloc(HASH_SIZE, uint16_image_color_table_node_t*); for (int32_t i = 0; i < HASH_SIZE; ++i) { cht[i] = NULL; } return cht; }