/* * This work is licensed under the Creative Commons Attribution 2.5 * Canada License. To view a copy of this license, visit * http://creativecommons.org/licenses/by/2.5/ca/ or send a letter to * Creative Commons, 171 Second Street, Suite 300, San Francisco, * California, 94105, USA. * * ********************************************************* * * UPSIZE: An image upsizing plugin featuring two modes. * * ***************************************************************************** * * 1. Smooth (IBFNBQH): Image upsizing with Interpolatory Box Filtered * * * Natural BiQuadratic Histosplines. * * * [This is the default option.] * * * * * * 2. Sharp (EANBQH): Exact Area image upsizing with Natural BiQuadratic * * * Histosplines. * * ***************************************************************************** * * For more details regarding these methods, see Fast Exact Area Image * Upsampling with Natural Biquadratic Histosplines by Nicolas * Robidoux, Adam Turcotte, Minglun Gong and Annie Tousignant, * pp.85-96 of Image Analysis and Recognition, 5th International * Conference, ICIAR 2008, Póvoa de Varzim, Portugal, June 25-27, * 2008. Proceedings, Aurelio C. Campilho, Mohamed S. Kamel (Eds.). * Lecture Notes in Computer Science 5112, Springer 2008, ISBN * 978-3-540-69811-1. Only the "Sharp" version is explicitly * discussed; "Smooth" differs in that it is box filtered (and in that * this implementation uses a different image size convention than the * "Sharp" resampler). */ /* * Upsize -- image resizing plug-in for The Gimp image manipulation * program Copyright (C) 2009 Nicolas Robidoux and Adam Turcotte * * LAST UPDATED: * March 25, 2011 * - Fixed Smooth/Sharp radio button callback. * - Uses gimp_image_insert_layer() for Gimp 2.7 and newer * instead of gimp_image_add_layer(), which has been deprecated. * * INSTALLATION INSTRUCTIONS: * Make sure the libgimp2.0-dev package is installed, then type: * gimptool-2.0 --install upsize.c * * NOTE: * This current version of the code has some unneeded * redundancies. Manana. */ /* * Natural boundary conditions mean that the normal derivative of the * histopolant (average matching reconstructed intensity profile) is * zero at the boundary. A alternate version of the code could use * not-a-knot boundary conditions instead (weaker variational * principle but would produce a scheme which is exact on affine * intensity profiles). */ /* * This source code can generate plug-ins which use either 16bit or * 8bit arrays to store the B-Spline coefficients after they have been * computed "in the vertical direction." * * More details: This tensor scheme requires one tridiagonal solve per * column, and one tridiagonal solve per row. The column solves are * performed first, with the results stored in integer arrays. The * results are retrieved when the row solves are performed. * * Storing the intermediate results in 8bit saves approximately half * of the memory compared to 16bit (not counting the input and output * images themselves), at the cost of introducing intensity errors * which are at most 1. * * If EIGHT_BIT_STORAGE is #defined, the half-computed B-spline * coefficients will be stored in an 8bit integer array. Otherwise, * they will be stored in a 16bit integer array. */ /* #define EIGHT_BIT_STORAGE */ /* * Increase these limits on the image dimensions as appropriate. To * our surprise, we have not encountered problems with index overflow * (in coarse_to_fine_coefficients, notably). */ #define MAX_WIDTH 65534 #define MAX_HEIGHT 65534 /* * This source code can generate executables which printf short * reports of time spent if the gimp was started from a terminal. * * If TIMING_REPORT is #defined, the time elapsed in the key part of * the program is sent to standard output. */ /* #define TIMING_REPORT */ /* * The ".5f +" is so that truncation (which the cast from float to int * effects for positive values) performs round to nearest. Negative * values are clamped later on to 0, so the fact that truncation * brings things "up" for negative results is not a concern. Halfway * points between nonnnegative integers are mapped to the larger * integer, which is acceptable---maybe even good if the image is * overexposed---but not important: round to nearest, with any * convention (banker's, toward zero---which may have advantages---or * toward +infty, which is the case now, would also be fine). Note * that an implicit cast to guchar is performed after clamping. */ #define F_ROUND_AND_CLAMP_0255(x) \ ({ \ gint rx = .5f + (x); \ guchar out = \ (rx>=0) ? ( (rx<=255) ? (guchar)rx : (guchar)255 ) : (guchar)0; \ out; \ }) /* * Key matrix factorization coefficients when computing B-Spline * coefficients: * * For natural boundary conditions: * C0 = 1/5 * C1 = 5/19 = 1 / ( 4 - C0 ) * C2 = 19/71 = 1 / ( 4 - C1 ) * C3 = 71/265 = 1 / ( 4 - C2 ) * C4 = 265/989 = 1 / ( 4 - C3 ) * C5 = 989/3691 = 1 / ( 4 - C4 ) * CINFTY = 2 - sqrt(3) = limit of this recurrence relation * CLAST = 1 / ( 5 - C_INFTY ) = ( 3 - sqrt(3) ) / 6 * * Convergence of the coefficients for natural boundary conditions: In * 32 bit single precision, C6 and all subsequent ones are equal to * CINFTY. Although C5 is different from CINFTY in this precision, it * differs from it only by 1.271371644 * 10^-7, just enough to be * different in standard single precision, but nonetheless an * insignificant amount (.2679491f instead of .2679492f). However, in * the interest of preventing integer overflow, we actually set use C5 * as a proxy for CINFTY in the computations which take place prior to * storing partically computed results into int8s. This is discussed * further below. */ #define C0_F .2000000f #define MINUS_C0_F -.2000000f #define C1_F .2631579f #define MINUS_C1_F -.2631579f #define C2_F .2676056f #define MINUS_C2_F -.2676056f #define C3_F .2679245f #define MINUS_C3_F -.2679245f #define C4_F .2679474f #define MINUS_C4_F -.2679474f #define C5_F .2679491f #define MINUS_C5_F -.2679491f #define CLAST_F .2113249f /* * In the interest of making sure that we don't overflow when * converting partially computed coefficients to integer data types, * in the computations which come prior to the cast we use C5 as a * proxy for CINFTY: C5 turns out to be the truncated value of CINFTY, * which is what we want. In the computations which take place after * the partially computed B-spline coefficients have been converted * back to floats, we use the correctly rounded value of CINFTY. */ #define CINFTY_F .2679492f #define MINUS_CINFTY_F -.2679492f /* * Implicitly, the code rescales the 0..255 uchar values so that they * fit in -255 to 255. Because the infinity norm of the inverse of the * tensor component of the "interpolation matrix" is 1/2, the computed * coefficients fit into -127.5 to 127.5, which we later encode into 0 * to 65535. * * The purpose of FIT_TO_RANGE is to turn pixel values ranging from 0 * to 255 to values ranging from -255 to 255. Values in this range are * mapped (arbitrarily close to when the matrix is large) values * strictly in the range -127.5 to 127.5 by the inverse of the matrix * A = * * 5 1 * 1 4 1 * 1 4 1 * 1 4 1 * 1 4 1 * ... * 1 5 * * which is the matrix for natural boundary conditions. * * FIT_TO_SCALE scales by a factor of 2 (and shifts). * * One can understand the factor of two as being part of solving when * the matrix is actually half the one shown above. This half-matrix * has an interpretation in terms of rescaling B-splines. * * The advantage of this interpretation is that this scaling * (asymptotically) keeps the range invariant: * * If all entries of y are within -127.5 to 127.5 (inclusive), and x * is the solution of A'x = y, then all entries of x are within the * same range (exclusive), provided * * A' = * 2.5 .5 * .5 2 .5 * .5 2 .5 * .5 2 .5 * ... * .5 2.5 * * Halving the matrix is equivalent to doubling the RHS. * * Yet a third interpretation: In order to store results using the * first version of the matrix, we double them (at any stage), and * then have to halve them when we want to use them. * * The important thing to keep in mind for later is that the values * are scaled by 2. */ #define FIT_TO_RANGE(x) ( 2 * ( (gint) (x) ) - 255 ) /* * The above can be rewritten using integer arithmetic, at the * cost of an additional (explicit) intermediate cast: */ /* #define FIT_TO_RANGE(x) ({gint in = (x); (in+in) - 255;}) */ #include #include #include #include /* Struct for scale parameters: */ typedef struct { gint width; gint height; gboolean mode; /* 0 = smooth, 1 = sharp */ } ScaleVals; /* Enum for reset button: */ enum { RESPONSE_RESET }; /* Use svals to access scale params. */ static ScaleVals svals; /* Dialog callback */ gdouble _upsize_width; gdouble _upsize_height; /* Aspect ratio */ gdouble _upsize_m_over_n; static void query (void); static void run (const gchar *name, gint nparams, const GimpParam *param, gint *nreturn_vals, GimpParam **return_vals); static gint32 scale_up_smooth (gint32 image_ID, GimpDrawable *drawable, gint m, gint n, gint mm, gint nn, gint channels); static void first_past_index (gint o, gint oo, gint first_past_kk[]); static void smooth_coarse_to_fine_coefficients (gint o, gint oo, gint first_past_kk[], gfloat left[], gfloat center[], gfloat right[], gfloat farright[]); static gint32 scale_up_sharp (gint32 image_ID, GimpDrawable *drawable, gint m, gint n, gint mm, gint nn, gint channels); static void last_overlapping_index (gint o, gint oo, gint last_overlapping_kk[]); static void coarse_to_fine_coefficients (gint o, gint oo, gint last_overlapping_kk[], gfloat left[], gfloat center[], gfloat right[], gfloat farright[]); static gboolean scale_dialog (gint32 image_ID, gint m, gint n); static void scale_callback (GtkWidget *widget, gpointer data); static void set_mode (GtkRadioButton *button, ScaleVals *svals ); GimpPlugInInfo PLUG_IN_INFO = { NULL, NULL, query, run }; MAIN() static void query (void) { gchar *help_path; gchar *help_uri; static GimpParamDef args[] = { { GIMP_PDB_INT32, "run-mode", "Run mode" }, { GIMP_PDB_IMAGE, "image", "Input image" }, { GIMP_PDB_DRAWABLE, "drawable", "Input drawable" }, { GIMP_PDB_INT32, "width", "Width of output" }, { GIMP_PDB_INT32, "height", "Height of output" } }; static GimpParamDef ret[] = { { GIMP_PDB_IMAGE, "image_out", "Output image" } }; /* use gimp_data_directory()/upsize-help for path */ help_path = g_build_filename (gimp_data_directory(), "upsize-help", NULL); help_uri = g_filename_to_uri (help_path, NULL, NULL); gimp_plugin_help_register (help_path, help_uri); g_free (help_path); g_free (help_uri); gimp_install_procedure ( "Upsize", "Upsize", "Image enlargement with biquadratic histosplines", "Nicolas Robidoux, Adam Turcotte", "Copyright Nicolas Robidoux, Adam Turcotte", "2007--2009", "_Upsize", "RGB*, GRAY*", GIMP_PLUGIN, G_N_ELEMENTS (args), G_N_ELEMENTS (ret), args, ret); gimp_plugin_menu_register ("Upsize", "/Image"); } static void run (const gchar *name, gint nparams, const GimpParam *param, gint *nreturn_vals, GimpParam **return_vals) { #ifdef TIMING_REPORT GTimer *timer; #endif static GimpParam values[2]; GimpPDBStatusType status = GIMP_PDB_SUCCESS; GimpRunMode run_mode; GimpDrawable *drawable; gint32 image_ID; gint32 new_image_ID; gint m; /* Height (in pixels) of the input image. */ gint n; /* Width (in pixels) of the input image. */ gint mm; /* Height (in pixels) of the output image. */ gint nn; /* Width (in pixels) of the output image. */ gint channels; /* Number of colour channels of the input image. */ /* * Set mandatory output values: */ *nreturn_vals = 2; *return_vals = values; values[0].type = GIMP_PDB_STATUS; values[0].data.d_status = status; /* * Getting run_mode (don't display a dialog in NONINTERACTIVE mode): */ run_mode = param[0].data.d_int32; image_ID = param[1].data.d_int32; /* * Get the specified drawable: */ drawable = gimp_drawable_get (param[2].data.d_drawable); _upsize_height = m = drawable->height; _upsize_width = n = drawable->width; channels = gimp_drawable_bpp (drawable->drawable_id); /* * Display warning if the input image is not large enough: */ if ( (m<7) || (n<8) ) { g_message ( "Only enlarges images at least\n8 pixels wide\nand 7 pixels tall.\n"); return; } switch (run_mode) { svals.mode = 1; case GIMP_RUN_INTERACTIVE: /* Initialize: */ svals.height = m; svals.width = n; _upsize_m_over_n = (gdouble) m / n; /* Display the dialog: */ if (! scale_dialog (image_ID, m, n)) return; mm = svals.height; nn = svals.width; /* * Set up a tile cache large enough to contain a column of input * tiles, or a row of output tiles: */ gimp_tile_cache_ntiles (((nn>m) ? nn : m) / gimp_tile_width () + 1); #ifdef TIMING_REPORT timer = g_timer_new (); #ifdef EIGHT_BIT_STORAGE g_print ("Upsize (with uint8 coefficient storage) is taking...\n"); #else g_print ("Upsize is taking...\n"); #endif #endif if (svals.mode) new_image_ID = scale_up_sharp (image_ID, drawable, m, n, mm, nn, channels); else new_image_ID = scale_up_smooth (image_ID, drawable, m, n, mm, nn, channels); #ifdef TIMING_REPORT g_print ("%g seconds to complete the computation.\n", g_timer_elapsed (timer, NULL)); g_timer_destroy (timer); #endif values[1].type = GIMP_PDB_IMAGE; values[1].data.d_image = new_image_ID; gimp_displays_flush (); /* flush to send data to core */ gimp_drawable_detach (drawable); /* then detach */ /* Set options in the core */ gimp_set_data ("Upsize", &svals, sizeof (ScaleVals)); if (new_image_ID != -1) gimp_display_new(new_image_ID); break; case GIMP_RUN_NONINTERACTIVE: if (nparams != 5) status = GIMP_PDB_CALLING_ERROR; if (status == GIMP_PDB_SUCCESS) { svals.width = param[3].data.d_int32; svals.height = param[4].data.d_int32; svals.mode = param[5].data.d_int32; mm = svals.height; nn = svals.width; /* * Set up a tile cache large enough to contain a column of input * tiles, or a row of output tiles: */ gimp_tile_cache_ntiles (((nn>m) ? nn : m) / gimp_tile_width () + 1); #ifdef TIMING_REPORT timer = g_timer_new (); #ifdef EIGHT_BIT_STORAGE g_print ("Upsize (with uint8 coefficient storage) is taking...\n"); #else g_print ("Upsize is taking...\n"); #endif #endif if (svals.mode) new_image_ID = scale_up_sharp (image_ID, drawable, m, n, mm, nn, channels); else new_image_ID = scale_up_smooth (image_ID, drawable, m, n, mm, nn, channels); #ifdef TIMING_REPORT g_print ("%g seconds to complete the computation.\n", g_timer_elapsed (timer, NULL)); g_timer_destroy (timer); #endif values[1].type = GIMP_PDB_IMAGE; values[1].data.d_image = new_image_ID; gimp_displays_flush (); /* flush to send data to core... */ gimp_drawable_detach (drawable); /* ...then detach */ /* Set options in the core */ gimp_set_data ("Upsize", &svals, sizeof (ScaleVals)); } break; case GIMP_RUN_WITH_LAST_VALS: gimp_get_data ("Upsize", &svals); mm = svals.height; nn = svals.width; /* * Set up a tile cache large enough to contain a column of input * tiles, or a row of output tiles: */ gimp_tile_cache_ntiles (((nn>m) ? nn : m) / gimp_tile_width () + 1); #ifdef TIMING_REPORT timer = g_timer_new (); #ifdef EIGHT_BIT_STORAGE g_print ("Upsize (with uint8 coefficient storage) is taking...\n"); #else g_print ("Upsize is taking...\n"); #endif #endif if (svals.mode) new_image_ID = scale_up_sharp (image_ID, drawable, m, n, mm, nn, channels); else new_image_ID = scale_up_smooth (image_ID, drawable, m, n, mm, nn, channels); #ifdef TIMING_REPORT g_print ("%g seconds to complete the computation.\n", g_timer_elapsed (timer, NULL)); g_timer_destroy (timer); #endif values[1].type = GIMP_PDB_IMAGE; values[1].data.d_image = new_image_ID; gimp_displays_flush (); /* flush to send data to core */ gimp_drawable_detach (drawable); /* then detach */ break; default: break; } return; } static gint32 scale_up_smooth (gint32 image_ID, GimpDrawable *drawable, gint m, gint n, gint mm, gint nn, gint channels) { gint first_past_ii[m-1]; gint first_past_jj[n-1]; gfloat top[mm]; gfloat middle[mm]; gfloat bottom[mm]; gfloat farbottom[mm]; gfloat left[nn]; gfloat center[nn]; gfloat right[nn]; gfloat farright[nn]; GimpDrawable *new_drawable; gint32 new_image_ID; gint32 new_layer_ID; GimpPixelRgn input_image, output_image; /* * Useful iterator bounds: */ gint nn_times_channels = nn * channels; gint n_times_channels = n * channels; gint n_minus_1_times_channels = (n-1) * channels; gint n_minus_7_times_channels = (n-7) * channels; gint n_minus_8_times_channels = n_minus_7_times_channels - channels; gint n_minus_2 = n - 2; gint m_times_channels = m * channels; gint m_minus_6_times_channels = (m - 6) * channels; gint m_minus_7_times_channels = m_minus_6_times_channels - channels; gint m_minus_2 = m - 2; /* * Utility iterators: */ gint k; /* Index related to the number of pixel values in an image row (or column). If the image has, say, three channels, to traverse a row one needs to visit 3n entries; k is used for that. */ gint c; /* Index (usually) related to the number of color channels (3 for RGB, 1 for Grayscale...). */ gint kk; /* * Row and column indices: */ gint i; /* Index of the input row under consideration. */ gint j; /* Index of the input column under consideration. */ gint ii; /* Index of the output row under consideration. */ gint last_ii; /* Local loop bound for ii. */ gint jj; /* Index of the output column under consideration. */ gint last_jj; /* Local loop bound for jj. */ /* * Computed coefficients/partial values: */ gfloat top_ii; gfloat middle_ii; gfloat bottom_ii; gfloat farbottom_ii; gfloat left_jj; gfloat center_jj; gfloat right_jj; gfloat farright_jj; #ifdef TIMING_REPORT GTimer *local_timer; gfloat local_time; #endif gfloat a_top_left[channels]; gfloat a_top_center[channels]; gfloat a_top_right[channels]; gfloat a_top_farright[channels]; gfloat a_middle_left[channels]; gfloat a_middle_center[channels]; gfloat a_middle_right[channels]; gfloat a_middle_farright[channels]; gfloat a_bottom_left[channels]; gfloat a_bottom_center[channels]; gfloat a_bottom_right[channels]; gfloat a_bottom_farright[channels]; gfloat a_farbottom_left[channels]; gfloat a_farbottom_center[channels]; gfloat a_farbottom_right[channels]; gfloat a_farbottom_farright[channels]; gfloat coef_left[channels]; gfloat coef_center[channels]; gfloat coef_right[channels]; gfloat coef_farright[channels]; /* * Computation of the B-spline coefficients: */ gfloat *a_row_p1, *a_row_p2; /* * Rescaled partially computed B-Spline coefficients are stored in * a. */ #ifdef EIGHT_BIT_STORAGE guchar *a, *a_ptr; a = g_new (guchar, n * m_times_channels); #else guint16 *a, *a_ptr; a = g_new (guint16, n * m_times_channels); #endif #ifdef TIMING_REPORT local_timer = g_timer_new (); local_time = g_timer_elapsed (local_timer, NULL); #endif first_past_index(m, mm, first_past_ii); first_past_index(n, nn, first_past_jj); #ifdef TIMING_REPORT g_print ("%g seconds for the first_past_index calls.\n", g_timer_elapsed (local_timer, NULL)-local_time); local_time = g_timer_elapsed (local_timer, NULL); #endif /* * Computation of the tensor components of the linear transformation * from B-spline coefficents to fine cell averages (actually fine * cell integrals, since we have already rescaled by the reciprocal * of the fine cell areas): */ smooth_coarse_to_fine_coefficients(m, mm, first_past_ii, top, middle, bottom, farbottom); smooth_coarse_to_fine_coefficients(n, nn, first_past_jj, left, center, right, farright); #ifdef TIMING_REPORT g_print ("%g seconds for the smooth_coarse_to_fine_coefficients calls.\n", g_timer_elapsed (local_timer, NULL)-local_time); #endif /* * Initialize the input pixel region: */ gimp_pixel_rgn_init (&input_image, drawable, 0, 0, n, m, FALSE, FALSE); /* * Set up an informative progress bar: */ gimp_progress_init ( "Upsizing with biquadratic histosplines"); gimp_progress_update(.005); #ifdef TIMING_REPORT local_time = g_timer_elapsed (local_timer, NULL); #endif { /* * Column-based (row by row) forward substitution, performed one * column at a time. */ guchar *input_col, *input_col_p; /* Storage and moving pointer for input image columns. */ gfloat *a_col, *a_col_p1, *a_col_p2; /* The computation is performed in the a_col array, accessed using the two corresponding pointers. */ input_col = g_new (guchar, m_times_channels); a_col = g_new (gfloat, m_times_channels); for (j=0; jdrawable_id), 100, GIMP_NORMAL_MODE); /* * Use gimp_image_insert_layer if version is at least 2.7.0. */ if (GIMP_MAJOR_VERSION >= 2 && GIMP_MINOR_VERSION >= 7) gimp_image_insert_layer (new_image_ID, new_layer_ID, -1, 0); else gimp_image_add_layer (new_image_ID, new_layer_ID, 0); new_drawable = gimp_drawable_get (new_layer_ID); /* * Init the output pixel region: */ gimp_pixel_rgn_init (&output_image, new_drawable, 0, 0, nn, mm, TRUE, TRUE); { gfloat *a_middle, *a_bottom, *a_farbottom, *a_top; gfloat *a_middle_p, *a_bottom_p, *a_farbottom_p, *a_top_p, *a_temp; /* * Storage for a pixel row to be output to new image buffer: */ guchar *output_row; guchar *output_row_ptr; gint nn_minus_1 = nn-1; a_top = g_new (gfloat, n_times_channels); a_middle = g_new (gfloat, n_times_channels); a_bottom = g_new (gfloat, n_times_channels); a_farbottom = g_new (gfloat, n_times_channels); /* * Temporary storage for a row of the output image: */ output_row = g_new (guchar, nn_times_channels); a_ptr = a; /* * ADAM AND NICOLAS: WHEN WE HAVE TIME, REWRITE SO NO FLOP IS * INVOLVED (using integer arithmetic): */ #ifdef EIGHT_BIT_STORAGE #define UNSCALING ( (gfloat) (.5) ) #define UNSCALED_SHIFT (-85) #else #define UNSCALING ( (gfloat) (1./514.) ) #define UNSCALED_SHIFT (-21845) #endif /* * Row solve on a_middle: */ a_row_p1 = a_middle; for (kk=n_times_channels; kk; kk--, a_row_p1++, a_ptr++) *a_row_p1 = ( (gint) *a_ptr + UNSCALED_SHIFT ) * UNSCALING; /* * Forward substitution: */ /* * We'll take care of the very first pixel on the way back. */ a_row_p1 = a_middle+channels; a_row_p2 = a_middle; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C0_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C1_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C2_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C3_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C4_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C5_F; /* * Because integer overflow is not a risk any more---we do clamp * values later on, and we'd have to even in exact arithmetic---we * use the rounded value of CINFTY instead of its truncated one. */ for (kk=n_minus_8_times_channels; kk; kk--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_CINFTY_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 = ( *a_row_p2 * MINUS_CINFTY_F + *a_row_p1 ) * CLAST_F; a_row_p1--; a_row_p2 = a_row_p1; for (c=channels; c; c--) a_row_p1--; for (kk=n_minus_7_times_channels; kk; kk--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * CINFTY_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C5_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C4_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C3_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C2_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C1_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C0_F; /* * Row solve on a_bottom: */ a_row_p1 = a_bottom; for (kk=n_times_channels; kk; kk--, a_row_p1++, a_ptr++) *a_row_p1 = ( (gint) *a_ptr + UNSCALED_SHIFT ) * UNSCALING; /* * Forward substitution: */ /* * We'll take care of the very first pixel on the way back. */ a_row_p1 = a_bottom+channels; a_row_p2 = a_bottom; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C0_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C1_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C2_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C3_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C4_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C5_F; /* * Because integer overflow is not a risk any more---we do clamp * values later on, and we'd have to even in exact arithmetic---we * use the rounded value of CINFTY instead of its truncated one. */ for (kk=n_minus_8_times_channels; kk; kk--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_CINFTY_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 = ( *a_row_p2 * MINUS_CINFTY_F + *a_row_p1 ) * CLAST_F; a_row_p1--; a_row_p2 = a_row_p1; for (c=channels; c; c--) a_row_p1--; for (kk=n_minus_7_times_channels; kk; kk--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * CINFTY_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C5_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C4_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C3_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C2_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C1_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C0_F; /* * Row solve on a_farbottom: */ a_row_p1 = a_farbottom; for (kk=n_times_channels; kk; kk--, a_row_p1++, a_ptr++) *a_row_p1 = ( (gint) *a_ptr + UNSCALED_SHIFT ) * UNSCALING; /* * Forward substitution: */ /* * We'll take care of the very first pixel on the way back. */ a_row_p1 = a_farbottom+channels; a_row_p2 = a_farbottom; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C0_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C1_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C2_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C3_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C4_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_C5_F; /* * Because integer overflow is not a risk any more---we do clamp * values later on, and we'd have to even in exact arithmetic---we * use the rounded value of CINFTY instead of its truncated one. */ for (kk=n_minus_8_times_channels; kk; kk--, a_row_p1++, a_row_p2++) *a_row_p1 += *a_row_p2 * MINUS_CINFTY_F; for (c=channels; c; c--, a_row_p1++, a_row_p2++) *a_row_p1 = ( *a_row_p2 * MINUS_CINFTY_F + *a_row_p1 ) * CLAST_F; a_row_p1--; a_row_p2 = a_row_p1; for (c=channels; c; c--) a_row_p1--; for (kk=n_minus_7_times_channels; kk; kk--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * CINFTY_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C5_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C4_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C3_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C2_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C1_F; for (c=channels; c; c--, a_row_p1--, a_row_p2--) *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C0_F; gimp_progress_update(.005+((gdouble)m*n)/((gdouble)mm*nn)*.15); last_ii = first_past_ii[0]; for (ii=0; iidrawable_id, FALSE); gimp_drawable_update (new_drawable->drawable_id, 0, 0, nn, mm); gimp_drawable_detach (new_drawable); return new_image_ID; } static void first_past_index (gint o, gint oo, gint first_past_kk[]) { /* * The comments use the variables which are relevant to the use of * this function in the horizontal direction. That is, n is o, nn is * oo, j is k, and jj is kk. * * This subprogram puts in first_past_jj[j], for j in 0..n-2, the * index of the first fine cell with a "center" which is not to the * left of the center of the j+1th coarse cell (the last cell has * index n-1, and so has no right neighbour). An alternate * description: first_past_jj[j] is the index of the first fine cell * with center >= j+3/2. This index matters for the following * reason: It is the index at which the rightmost relevant B-Spline * coefficient shifts to the right (center = j+3/2 could go either * way). * * This code uses a general convention that the input pixels have * unit sides, and that (describing things in 1D) the global * coordinate of the left endpoint of an input cell is exactly its * index. That is, the first coarse (input) cell goes from 0 to 1, * the last coarse cell from n-1 to n. * * Because of this convention, there is an alternate description of * first_past_jj[j]: Because coarse cell centers are 1/2 to the * right of left boundaries, first_past_jj is also the index of the * first fine cell for which the left end of the box filtering mask * (which has diameter 1/2) is >= j+1. The very first "box" starts * at zero, and it shifts by h at every index shift. This gives jj*h * >= j+1, that is, it is the first jj which violates jj*(n-1) < * (j+1)*(nn-1). * * Because the jth cell is [j,j+1], and choosing the very first fine * cell center to be at 1/2 (same as the very first coarse cell * center), so that the very last fine, and coarse, cell center is * n-1/2, we have that as the fine index runs from 0 to nn-1, a * distance of n-1 is covered, so that * * h = (n-1)/(nn-1) * * (the usual step for interpolation methods). Consequently, * first_past_jj[j] is the first index jj for which * * 1/2 + jj * h >= j + 3/2; * * in integer arithmetic: * * jj * (n-1) >= (j+1) * (nn-1). * * Alternately, it is the first index for which the condition * * jj * (n-1) < (j+1) * (nn-1) * * does not hold. */ gint o_minus_1 = o-1; gint k = 0; if (oo>o) { gint kk = 1; gint kk_times_o_minus_1 = o_minus_1; gint oo_minus_1 = oo-1; gint k_plus_1_times_oo_minus_1 = oo_minus_1; do { while (kk_times_o_minus_1 < k_plus_1_times_oo_minus_1) { kk++; kk_times_o_minus_1 += o_minus_1; } first_past_kk[k]=kk; kk++; kk_times_o_minus_1 += o_minus_1; k++; k_plus_1_times_oo_minus_1 += oo_minus_1; } while (kdrawable_id), 100, GIMP_NORMAL_MODE); /* * Use gimp_image_insert_layer if version is at least 2.7.0. */ if (GIMP_MAJOR_VERSION >= 2 && GIMP_MINOR_VERSION >= 7) gimp_image_insert_layer (new_image_ID, new_layer_ID, -1, 0); else gimp_image_add_layer (new_image_ID, new_layer_ID, 0); new_drawable = gimp_drawable_get (new_layer_ID); /* * Init the output pixel region: */ gimp_pixel_rgn_init (&output_image, new_drawable, 0, 0, nn, mm, TRUE, TRUE); /* * Given that we will not be storing results back into guchars * until the very end, we start by undoing the effect of * shifting by .5 to store -32767.5..32767.5 in int16, as well * as of the FIT_TO_RANGE rescaling. (If this was in 3D, we * would postpone this "unscaling" until dealing with the very * last direction, or we would allocate a float "slab" instead * of a row.) * * Once this is done, we can compute solutions taking results at * face value. * * We also divide by the fine cell areas, which allows us to * compute with integrals instead of averages later on. * * NATURAL BOUNDARY CONDITIONS: * * Given the solution z of * * A ( m2 * z + b2 ) = m1 * y + b1, * * from which one gets that * * y = A ( m2/m1 * z + b2/m1 ) - b1/m1, * * so that * * A^{-1} y = m2/m1 * z + b2/m1 - b1/m1 * A^{-1} 1 * * This gives us the solution of A x = y. * * For natural boundary conditions, m1 = 2, b1 = -255. * * Now, z is the value stored as a uint16, so that m2 * z + b2 * represents the value prior the conversion. Since the conversion * is * * z = 257 * a + 32767 + 1/2 * * we have that the solution a we had prior to normalization was * * ( 2z - 65535 )/514, so that m2 = 1/257 and b2 = -255/2. * * It is easy to see that A^{-1} 1 = 1/6 when A = * * 5 1 * 1 4 1 * 1 4 1 * 1 4 1 * 1 4 1 * ... * 1 5 * * (Here, m is a scalar, and b and c are constant shifts. There * is some abuse of notation in the above.) * * This implies that the shift if -85/2 (=b2/m1 - b1/m1 * 1/6) * and the scaling is m2/m1 = 1/514. * * If we are using uchars for storage, things are a bit * different: * * m1 and b1 are as before, but m2 = 1 and b2 = -127.5. This * gives shift = -85/2 and scaling = 1/2. * * The corresponding constants, unscaling and unscaled_shift are * defined in the calling program just prior to the first call. */ a_ptr = a; /* * Row solve on a_middle: */ a_row_p1 = a_middle; c = n_times_channels; do { *a_row_p1++ = ( (gint) *a_ptr++ + UNSCALED_SHIFT ) * unscaling; } while (--c); /* * Forward substitution: */ /* * We'll take care of the very first pixel on the way back. */ a_row_p2 = a_middle; a_row_p1 = a_middle+channels; c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C0_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C1_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C2_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C3_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C4_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C5_F; } while (--c); /* * Because integer overflow is not a concern any more---we do * clamp values later on, and we'd have to even in exact * arithmetic---we use the rounded value of CINFTY instead of * its truncated one. */ c = n_minus_8_times_channels; while (c--) { *a_row_p1++ += *a_row_p2++ * MINUS_CINFTY_F; } c = channels; do { *a_row_p1 = ( *a_row_p2++ * MINUS_CINFTY_F + *a_row_p1 ) * CLAST_F; ++a_row_p1; } while (--c); a_row_p2 = --a_row_p1; c = channels; do { --a_row_p1; } while (--c); c = n_minus_7_times_channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * CINFTY_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C5_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C4_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C3_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C2_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C1_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C0_F; --a_row_p1; } while (--c); /* * Row solve on a_bottom: */ a_row_p1 = a_bottom; c = n_times_channels; do { *a_row_p1++ = ( (gint) *a_ptr++ + UNSCALED_SHIFT ) * unscaling; } while (--c); a_row_p2 = a_bottom; a_row_p1 = a_bottom+channels; c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C0_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C1_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C2_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C3_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C4_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C5_F; } while (--c); c = n_minus_8_times_channels; while (c--) { *a_row_p1++ += *a_row_p2++ * MINUS_CINFTY_F; } c = channels; do { *a_row_p1 = ( *a_row_p2++ * MINUS_CINFTY_F + *a_row_p1 ) * CLAST_F; ++a_row_p1; } while (--c); a_row_p2 = --a_row_p1; c = channels; do { --a_row_p1; } while (--c); c = n_minus_7_times_channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * CINFTY_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C5_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C4_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C3_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C2_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C1_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C0_F; --a_row_p1; } while (--c); /* * Row solve on a_farbottom: */ a_row_p1 = a_farbottom; c = n_times_channels; do { *a_row_p1++ = ( (gint) *a_ptr++ + UNSCALED_SHIFT ) * unscaling; } while (--c); a_row_p2 = a_farbottom; a_row_p1 = a_farbottom+channels; c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C0_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C1_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C2_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C3_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C4_F; } while (--c); c = channels; do { *a_row_p1++ += *a_row_p2++ * MINUS_C5_F; } while (--c); c = n_minus_8_times_channels; while (c--) { *a_row_p1++ += *a_row_p2++ * MINUS_CINFTY_F; } c = channels; do { *a_row_p1 = ( *a_row_p2++ * MINUS_CINFTY_F + *a_row_p1 ) * CLAST_F; ++a_row_p1; } while (--c); a_row_p2 = --a_row_p1; c = channels; do { --a_row_p1; } while (--c); c = n_minus_7_times_channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * CINFTY_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C5_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C4_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C3_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C2_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C1_F; --a_row_p1; } while (--c); c = channels; do { *a_row_p1 = ( *a_row_p1 - *a_row_p2-- ) * C0_F; --a_row_p1; } while (--c); last_ii = last_overlapping_ii[0]; for (ii=0; iidrawable_id, FALSE); gimp_drawable_update (new_drawable->drawable_id, 0, 0, nn, mm); gimp_drawable_detach (new_drawable); return new_image_ID; } static void last_overlapping_index (gint o, gint oo, gint last_overlapping_kk[]) { /* * The comments use the variables which are relevant to the use of * this function in the horizontal direction. That is, n is o, nn is * oo, j is k, and jj is kk. * * This code uses a general convention that the input pixels have * unit sides, and that (describing things in 1D) the global * coordinate of the left endpoint of an input cell is exactly its * index. That, is the first input cell goes from 0 to 1, the last * one from n-1 to n. Given that, * * jj * n / nn (in rational arithmetic) * * is the absolute coordinate of the left endpoint of the index jj * fine cell (recall that jj starts at zero). Given that the left * endpoints j of the coarse cells are 0, 1, 2... this implies that * * (j * nn ) / n (in integer arithmetic) * * is the index of the first fine cell which is contained in the * index j coarse cell, so that * * ((j+1) * nn ) / n - 1 (in integer arithmetic) * * is the index of the last fine cell which overlaps the index j * coarse cell. * * There is another way of caracterizing this index: It is the * largest jj such that jj * n < j + 1. Using this caracterization * allows one to avoid integer divisions, and because by sweeping * through the jjs we know approximately what this index is, few * iterations are necessary. The key idea is as follows: whenever * the product of the next candidate jj by n initially fails the * above condition, the current jj is last_jj for this value of j. * * If n=nn, the last overlapping index is equal to the index. */ gint o_minus_1 = o-1; /* Useful iteration bound. */ gint k = 0; if (oo>o) { gint k_plus_one_times_oo = oo; gint kk = 0; gint kk_plus_one_times_o = o; do { ++kk; /* Because 1 < n <= nn, we know that the first overlapping jj can't be the last overlapping jj, so we can increment without checking the condition. */ kk_plus_one_times_o += o; while (kk_plus_one_times_o < k_plus_one_times_oo) { ++kk; /* If the next jj satisfies the key condition, make it jj. */ kk_plus_one_times_o += o; } last_overlapping_kk[k] = kk; /* We have found the last overlapping jj for the current value of j. */ k_plus_one_times_oo += oo; } while (++k