/* Last edited on 2024-08-31 06:02:32 by stolfi */ /* ---------------------------------------------------------------------- */ #define PROG_HELP := \ " [ -inverse ] \\\n" \ #define PROG_INFO := \ " -inverse\n" \ " If present, this option requests the inverse of the" \ " projective map indicated by the \"-matrix\" or" \ " \"-points\" options. (For the latter, it should" \ " be equivalent to swapping the to wpoint lists.).\n" \ "\n" \ "MODIFICATION HISTORY\n" \ " jan/2008 Added the \"-inverse\" option. J. Stolfi, IC-UNICAMP.\n" \ typedef struct options_t { bool_t inverse; /* Requests the inverse map. */ } options_t; void get_options(void) { o->inverse = argparser_keyword_present(pp, "-inverse"); } /* ---------------------------------------------------------------------- */ void process_map_unmap_commands ( char *mapFile, mpoint_vec_t *mpoints, hr2_pmap_t *isys, hr2_pmap_t *M, hr2_pmap_t *osys ); /* Processes the \"-map\" and \"-unmap\" options specified in {mpoints}, writing the result on the file named "{mapFile}", or to {stderr} if {mapFile} is "-". Uses the projective maps {isys} (pixel to user input coords), {M} (input to output user coords), and {osys} (pixel to user output coords). */ void process_map_unmap_commands ( char *mapFile, mpoint_vec_t *mpoints, hr2_pmap_t *isys, hr2_pmap_t *M, hr2_pmap_t *osys ) { FILE *wr = (strcmp(mapFile, "-") == 0 ? stderr : open_write(mapFile, TRUE)); for (int32_t k = 0; k < mpoints->ne; k++) { mpoint_t *mpk = &(mpoints->e[k]); r2_t u_in, u_out; /* User coordinates in the two images. */ r2_t p_in, p_out; /* Pixel coordinates in the two images. */ if (mpk->inv) { u_out = mpk->p; r2_map_projective(&u_out, &(M->inv), &u_in, NULL); } else { u_in = mpk->p; r2_map_projective(&u_in, &(M->dir), &u_out, NULL); } r2_map_projective(&u_in, &(isys->inv), &p_in, NULL); r2_map_projective(&u_out, &(osys->inv), &p_out, NULL); fprintf(wr, "%-8s", mpk->tag); fprintf(wr, " %7.1f %7.1f", u_in.c[0], u_in.c[1]); fprintf(wr, " %7.1f %7.1f", p_in.c[0], p_in.c[1]); fprintf(wr, " %7.1f %7.1f", p_out.c[0], p_out.c[1]); fprintf(wr, " %7.1f %7.1f", u_out.c[0], u_out.c[1]); fprintf(wr, "\n"); } fclose(wr); } /* Map/unmap given points: */ process_map_unmap_commands(o->mapFile, &(o->mpoints), &isys, &(o->M), &osys); void hr2_pmap_opt_decode_map(int32_t nv, double y[], hr2_pmap_t *M, double rad[]) { switch(type) { case hr2_pmat_type_TRANSLATION: { assert(nv == 2); /* The two variables are the displacement, with hyperbolic */ /* contraction so that the unit disk of {y} represents the disk */ /* of {\RR^2} with radius {r_max}. Undo the contraction: */ r2_t p = (r2_t){{ y[0], y[1] }}; r2_hyper_expand(&p, r_max, &p); (*M) = hr2_pmap_translation(&p); break; } case hr2_pmat_type_CONGRUENCE: { assert(nv == 3); /* The first two {y} variables are the displacement of the origin, as in translation. The thurd one is the angle mapped from [-2\pi _ +2\pi] to {[-1 _ +1]}. */ demand(FALSE< "** CONGRUENCE map not implemented"); break; } case hr2_pmat_type_SIMILARITY: { assert(nv == 4); /* The first two {y} variables are the displacement of the origin, as in translation. The other two are the image of the {(1,0)} vector, also hyper contracted. */ r2_t p = (r2_t){{ y[0], y[1] }}; r2_hyper_expand(&p, r_max, &p); r2_t u = (r2_t){{ y[2], y[3] }}; r2_hyper_expand(&u, r_max, &u); r2_t q; r2_add(&p, &u, &q); (*M) = hr2_pmap_similarity_from_two_points(&p, &q, flip); break; } case hr2_pmat_type_AFFINE: { assert(nv == 6); /* The first two {y} variables are the displacement of the origin, as in translation. The other four are the images of the {(1,0)} and {(0,1)} vectors, also hyper contracted. */ r2_t p = (r2_t){{ y[0], y[1] }}; r2_hyper_expand(&p, r_max, &p); r2_t u = (r2_t){{ y[2], y[3] }}; r2_hyper_expand(&u, r_max, &u); r2_t v = (r2_t){{ y[4], y[5] }}; r2_hyper_expand(&v, r_max, &v); r2x2_t E = (r2x2_t){{{ u[0], u[1] }, { v[0], v[1] }}}; (*M) = hr2_pmap_aff_from_mat_and_disp(&E, &p); break; } case hr2_pmat_type_PROJECTIVE: { assert(nv == 8); /* The point {h[0..7]} of {\RR^8} is scaled by {r_max} then mapped to the matrix {M.dir} viewed as a point of the units sphere in {\RR^9}, by stereographic projection from the matrix {M0} that has {M[0][0] = 1} and all other elements equal to zero; that is, scaling of {\RR^2} by {(0,0)}. Thus if {y} is inside {M} will never get too close to {M0}. */ double d2 = r_max*r*max*rn_norm_sqr(nv, y); break; } default: demand(FALSE, "unimplemented map type") } } /* ---------------------------------------------------------------------- */ // if (argparser_keyword_present(pp, "-ifilter")) // { o->nfilter_in = argparser_get_next_int(pp, 1, MAX_NFILTER); } // else // { o->nfilter_in = 2; } /* ---------------------------------------------------------------------- */ r2_t map_user_pixel(double col, double row, hr2_pmap_t *csys); /* Given the user indices {col,row} of a pixel, and the native-to-user projective map {csys}, returns the native coordinates of the pixel's center. */ r2_t map_user_pixel(double col, double row, hr2_pmap_t *csys) { /* Convert debug pixel indices to homogeneous coordinates {hp}: */ r3_t hp = (r3_t) {{ 1.0, col + 0.5, row + 0.5 }}; /* Map {hp} to {hq}, get the homogeneous coordinates {qw,qx,qy}: */ r3_t hq; r3x3_map_row(&hp, &(csys->inv), &hq); double qw = hq.c[0]; double qx = hq.c[1]; double qy = hq.c[2]; /* Check for validity of result: */ if (qw <= 0.0) { /* Result is at infinity, or beyond: */ return (r2_t){{ -1, -1 }}; } else { /* Convert {qw,qx,qy} to Cartesian coordinates {qX,qY}: */ double qX = qx/qw; double qY = qy/qw; return (r2_t){{ qX, qY }}; } } /* ---------------------------------------------------------------------- */ { /* Check for identity: */ T->ident_proj = hr2_pmap_is_ident(&(T->pmap)); o->T.ident_proj = TRUE; if (coef != el[0]*((double)(i==j))) { o->T.ident_proj = FALSE; } /* Compute matrix of inverse map: */ if (o->T.ident_proj) { o->T.PINV = o->T.PDIR; } else { } o->T.ident_proj = TRUE; /* Compute matrix of inverse map: */ if (o->T.ident_proj) { o->T.PINV = o->T.PDIR; } else { } } o->T.ident_proj = TRUE; } void fitr_transform_point ( double xi, double yi, fitr_transform_t *T, bool_t inv, double *woP, double *xoP, double *yoP, double *JP, bool_t debug ); /* Applies the transformation {T} to the pont {p=(xi,yi)}, returns the homogeneous coordinates in {*woP}, {*xoP}, {*yoP}. Also computes the Jacobian {*JP} of {T} at that point. If {inv}, uses the inverse of the map {T}. */ void fitr_apply_radial_map(double *wP, double *xP, double *yP, double R, double *J); /* Applies the radial map with parameter R to the point with homogeneous coordinates {[*wP, *xP, *yP]}. Also post-multiplies the 2x2 matrix {*J} by the Jacobian of the map. */ void fitr_apply_projective_map(fitr_projmap_t *M, double *wP, double *xP, double *yP, double *J); /* Applies the projective map {M} to the signed homogeneous vector {p = [*wP, *xP, *yP]}. The vector {[*wP, *xP, *yP]} is then normalized so that {*wP} is in {{-1,0,+1}}. Also post-multiplies the 2x2 matrix {*J} by the Jacobian of the map at {p}. */ " [ -unit { pixel | image } ] \\\n" \ "In both systems, coordinates are measured" \ " either in pixels, or in multiples of the relevant image" \ " dimension, depending on the \"-unit\" argument." " -unit { pixel | image }\n" \ " This optional argument specifies the unit of measurement" \ " for the input and output coordinate systems in both" \ " images: whether pixel widths and heights (\"pixel\"), or" \ " the image dimensions along the corresponding axis (\"image\"). With" \ " the second option, the image is assumed to have unit width and" \ " height, independently of its size in pixels. If omitted," \ " the program assumes \"-unit pixel\".\n" \ "\n" \ "\n" \ typedef struct fitr_projmap_t { /* Note: the order of fields is important for procedures below! */ double WT, XT, YT, WX, XX, YX, WY, XY, YY; } fitr_projmap_t; /* Coefficients of the planar projetive transformation `(X,Y) => (X',Y')' defined by the formulas W' = (WT + WX*X + WY*Y), X' = (XT + XX*X + XY*Y)/W', Y' = (YT + YX*X + YY*Y)/W', Note that the origin (0,0) is taken to (XT/WT, YT/WT), the infinity point of the X axis (+oo,0) goes to (XX/WX, YX/WX), the infinity point on the Y axis (0,+oo)goes to (XY/WY, YY/WY), and the local deformation (Jacobian) matrix near the origin is dX'/dX = XX/WT, dX'/dY = XY/WT, dY'/dX = YX/WT, dY'/dY = YY/WT. */ typedef struct transform_t { double ocol_in; /* Origin X pos in input image. */ double orow_in; /* Origin Y pos in input image. */ double R; /* Radial deformation parameter (0 = none). */ fitr_projmap_t PDIR; /* Projective map from input domain to output domain. */ fitr_projmap_t PINV; /* Inverse of {PDIR} */ bool_t ident_proj; /* TRUE if {PDIR} and {PINV} are the identity */ double ocol_ot; /* Origin X pos in output image. */ double orow_ot; /* Origin Y pos in output image. */ } transform_t; /* The transformation assumes that each pixel is a unit square with integer corners, and also that the origin of the input coordinate system lies `ocol_in' pixels to the right of the left edge, and `orow_in' pixels below the top edge of the input image. The fields `ocol_ot' and `orow_ot' similarly define the origin position in the output image. The transformation to be applied to the domain is a radial correction with parameter `T.R' relative to the origin, followed by a projective map with homogeneous matrix `T.P'. The latter assumes that the Y axis points UP, even though the parameters `orow_in' and `orow_ot' are measured from the top edge down. Thus the sample with indices `[row][col]' is assumed to be centered at point `((col + 0.5) - ocol, orow - (row + 0.5))'. If `T.R > 0', the radial transform scales each point at distance `d' by the factor `1/(1 - s2)' where `s2 = (d/T.R)^2'. For small d this is essentially `1 + s2'. As `d' approaches `T.R', the factor goes to +oo. If `T.R < 0', the inverse map is used; namely, each point at distance `d' is scaled by the factor `2/(1 + sqrt(1 + 4*s2))'. */ /* Provide default scale for input pixels: */ if (zero_in == +INF) { zero_in = 0; } if (unit_in == +INF) { unit_in = (double)maxval_in; } if (maxval_us == +INF) { maxval_us = (double)maxval_in; } hd_ot.forceplain = FALSE; switch (hd_in.format) { case PPM_FORMAT: case RPPM_FORMAT: if (maxval_us > PPM_MAXMAXVAL) { pm_error("maxval must be at most %d", PPM_MAXMAXVAL); } hd_ot.maxval = (int)maxval_us; hd_ot.format = hd_in.format; break; case PGM_FORMAT: case RPGM_FORMAT: if (maxval_us > PGM_OVERALLMAXVAL) { pm_error("maxval must be at most %d", PGM_OVERALLMAXVAL); } hd_ot.format = hd_in.format; hd_ot.maxval = (int)maxval_us; break; case PBM_FORMAT: case RPBM_FORMAT: if (maxval_us > PGM_OVERALLMAXVAL) { pm_error("maxval must be at most %d", PGM_OVERALLMAXVAL); } hd_ot.maxval = (int)maxval_us; if (hd_ot.maxval > 1) { pm_message("promoting input file to PGM"); /* The new pbmplus may write raw PGM files even when maxval > 255. */ /* Such files are incompatible with older versions of Netpbm; avoid them. */ if ((hd_in.format == PBM_FORMAT) || (hd_ot.maxval > PGM_MAXMAXVAL)) { hd_ot.format = PGM_FORMAT; hd_ot.forceplain = TRUE; } else { hd_ot.format = RPGM_FORMAT; } } break; default: pm_error("bad input file format"); } /* Provide default scale for output pixels: */ if (zero_ot == +INF) { zero_ot = 0; } if (unit_ot == +INF) { unit_ot = (double)hd_ot.maxval; } void compute_output_dims(int cols_in, int rows_in, double Xc, double Yc, double Rd, ProjMap *PP, int *cols_otP, int *rows_otP) /* Computes the output image size `*cols_otP' and `*rows_otP' so as to include all mapped parts of the input image, except possibly for parts that are mapped to negative X (off the left edge) or negative Y (off the bottom edge). /* { /* Compute coordinate ranges relative to the center point (Xc,Yc): */ double xmin = -Xc; double xmax = (double)cols_in - Xc; double ymin = - Yc; double ymax = (double)rows_in - Yc; /* Expand ranges so as to account for radial distortion: */ { double xp, yp, xq, yq; double xqmin = +INF, xqmax = -INF; double yqmin = +INF, yqmax = -INF; int ix, iy; /* Assumes that the radial distortion is monotonic in each quadrant. */ /* Thus the H-V bounding box is defined by the corners and axial pts. */ for (ix = -1; ix <= +1; ix++) for (iy = -1; iy <= +1; iy++) { double xp = (ix == 0 ? 0 : (ix < 0 ? xmin : xmax)); double yp = (iy == 0 ? 0 : (iy < 0 ? ymin : ymax)); if ( ((ix != 0) || (iy != 0)) && (xp >= xmin) && (xp <= xmax) && (yp >= ymin) && (yp <= ymax)) { } /* ---------------------------------------------------------------------- */ " -scale {SCALE}\n" \ " If this parameter is present, it applies the scaling" \ " factor {SCALE} to the output image, including the the" \ " size of the output image, the pixel colum and row of" \ " its coordinate system origin, and the projective" \ " map. If omitted, no extra scaling is applied (same as \"-scale 1\").\n" \ "\n" \ if (argparser_keyword_present(pp, "-scale")) { o->scale = argparser_get_next_double(pp, 0.001, 100.0); } else { o->scale = 1; } /* ---------------------------------------------------------------------- */ /* Map center of debugging pixel from pixel coords to pixel output coords {dbp}: */ r2_t dbp; bool_t dbp_defined; if ((o->debug.c[0] < 0) || (o->debug.c[1] < 0)) { dbp = (r2_t){{ NAN, NAN }}; dbp_defined = FALSE; } else { dbp.c[0] = o->debug.c[0]; dbp.c[1] = oRows - o->debug.c[1]; r2_t odbp; r2_map_projective(&dbp, &(osys.dir), &odbp, NULL); dbp_defined = TRUE; if (o->verbose) { fprintf(stderr, "watching output point with coordinates:\n"); fprintf(stderr, " (%7d,%7d) PBM (Y down)\n", o->debug.c[0], o->debug.c[1]); fprintf(stderr, " (%7.1f,%7.1f) libimg (Y up)\n", dbp.c[0], dbp.c[1]); fprintf(stderr, " (%7.1f,%7.1f) user\n", odbp.c[0], odbp.c[1]); } } void extract_matching_features ( feature_vec_t *fta, int32_vec_t *ixa, feature_vec_t *ftb, int32_vec_t *ixb, r2_t *fpa[], r2_t *fpb[] ); /* Assumes that {ixa} and {ixb} are vectors of the same length {nf} such that features {fta.e[ixa.e[kf]]} and {ftb.e[ixb.e[kf]]} have the same tag, for {kf} in {0..nf-1}. Stores into {fpa[0..nf-1]} the addresses of {fta.e[ixa.e[0..nf-1]].pos} and into {fpb[0..nf-1]} the addresses of {ftb.e[ixb.e[0..nf-1]].pos}. */ void extract_matching_features ( feature_vec_t *fta, int32_vec_t *ixa, feature_vec_t *ftb, int32_vec_t *ixb, r2_t *fpa[], r2_t *fpb[] ) { demand(ixa->ne == ixb->ne, "mismatched lengths {ixa,ixb}"); /* Paranoia. */ int32_t nf = ixa->ne; for (int32_t kf = 0; kf < nf; kf++) { fpa[kf] = &(fta->e[ixa->e[kf]].pos); fpb[kf] = &(ftb->e[ixb->e[kf]].pos); } } }