C 16
IBFNBQH.c Guest on 27th September 2020 05:40:17 PM
  1. /*
  2.  * IBFNBQH: Average matching image upsizing with Interpolatory Box Filtered
  3.  *          Natural Biquadratic Histosplines.
  4.  */
  5.  
  6. /* GIMP - The GNU Image Manipulation Program
  7.  * Copyright (C) 1995 Spencer Kimball and Peter Mattis
  8.  *
  9.  * IBFNBQH -- image resizing plug-in for The Gimp image manipulation
  10.  * program Copyright (C) 2007 Nicolas Robidoux and Adam Turcotte
  11.  *
  12.  * This program is free software; you can redistribute it and/or
  13.  * modify it under the terms of the GNU General Public License as
  14.  * published by the Free Software Foundation; either version 2 of the
  15.  * License, or (at your option) any later version.
  16.  *
  17.  * This program is distributed in the hope that it will be useful, but
  18.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  19.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20.  * General Public License for more details.
  21.  *
  22.  * You should have received a copy of the GNU General Public License
  23.  * along with this program; if not, write to the Free Software
  24.  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  25.  * 02111-1307, USA.
  26.  */
  27.  
  28. /*
  29.  * Natural boundary conditions mean that the normal derivative of the
  30.  * histopolant (average matching reconstructed intensity profile) is
  31.  * zero at the boundary. A future version of the code may use
  32.  * not-a-knot boundary conditions instead (weaker variational
  33.  * principle but would produce a scheme which is exact on affine
  34.  * intentisity profiles).
  35.  */
  36.  
  37. /*
  38.  * This source code can generate plug-ins which use either 16bit or
  39.  * 8bit arrays to store the B-Spline coefficients after they have been
  40.  * computed "in the vertical direction."
  41.  *
  42.  * More details: This tensor scheme requires one tridiagonal solve per
  43.  * column, and one tridiagonal solve per row. The column solves are
  44.  * performed first, with the results stored in integer arrays. The
  45.  * results are retrieved when the row solves are performed.
  46.  *
  47.  * Storing the intermediate results in 8bit saves approximately half
  48.  * of the memory compared to 16bit (not counting the input and output
  49.  * images themselves), at the cost of introducing intensity errors
  50.  * which are at most 1.
  51.  */
  52.  
  53. /*
  54.  * If EIGHT_BIT_STORAGE is #defined, the half-computed B-spline
  55.  * coefficients will be stored in an 8bit integer array. Otherwise,
  56.  * they will be stored in a 16bit integer array.
  57.  */
  58.  
  59. /* #define EIGHT_BIT_STORAGE */
  60.  
  61. /*
  62.  * This source code can generate executables which printf short
  63.  * reports of computation time if the gimp is started from a terminal.
  64.  *
  65.  * If TIMING_REPORT is #defined, these reports are created and sent to
  66.  * standard output.
  67.  */
  68.  
  69. /* #define TIMING_REPORT */
  70.  
  71. /*
  72.  * Increase these limits on the image dimensions as appropriate. To
  73.  * our surprise, we have not encountered problems with index overflow
  74.  * (in smooth_coarse_to_fine_coefficients, notably).
  75.  */
  76. #define MAX_WIDTH  65534
  77. #define MAX_HEIGHT 65534
  78.  
  79. #include <libgimp/gimp.h>
  80. #include <libgimp/gimpui.h>
  81. #include <math.h>
  82. #include <stdlib.h>
  83.  
  84. /* Enum for reset button */
  85. enum {
  86.     RESPONSE_RESET
  87. };
  88.  
  89. /* Struct for scale parameters: */
  90. typedef struct
  91. {
  92.   gint width;
  93.   gint height;
  94. } ScaleVals;
  95.  
  96. /* Dialog callback */
  97. gdouble _ibfnbqh_width;
  98. gdouble _ibfnbqh_height;
  99.  
  100. /* Aspect ratio */
  101. gdouble _ibfnbqh_m_over_n;
  102.  
  103. /* Use svals to access scale params. */
  104. static ScaleVals svals;
  105.  
  106. static void query (void);
  107.  
  108. static void run (const gchar      *name,
  109.                  gint              nparams,
  110.                  const GimpParam  *param,
  111.                  gint             *nreturn_vals,
  112.                  GimpParam       **return_vals);
  113.  
  114. static gint32 scale_up (gint32        image_ID,
  115.                         GimpDrawable *drawable,
  116.                         gint          m,
  117.                         gint          n,
  118.                         gint          mm,
  119.                         gint          nn,
  120.                         gint          channels);
  121.  
  122. static void first_past_index (gint o,
  123.                               gint oo,
  124.                               gint first_past_kk[]);
  125.  
  126. static void smooth_coarse_to_fine_coefficients (gint   o,
  127.                                                 gint   oo,
  128.                                                 gint   first_past_kk[],
  129.                                                 gfloat left[],
  130.                                                 gfloat center[],
  131.                                                 gfloat right[],
  132.                                                 gfloat farright[]);
  133.  
  134. static gboolean scale_dialog (gint32 image_ID,
  135.                               gint   m,
  136.                               gint   n);
  137.  
  138. static void scale_callback (GtkWidget *widget,
  139.                             gpointer   data);
  140.  
  141. GimpPlugInInfo PLUG_IN_INFO =
  142. {
  143.   NULL,
  144.   NULL,
  145.   query,
  146.   run
  147. };
  148.  
  149. MAIN()
  150.  
  151. static void
  152. query (void)
  153. {
  154.   gchar *help_path;
  155.   gchar *help_uri;
  156.  
  157.   static GimpParamDef args[] =
  158.   {
  159.     { GIMP_PDB_INT32,    "run-mode",  "Run mode"         },
  160.     { GIMP_PDB_IMAGE,    "image",     "Input image"      },
  161.     { GIMP_PDB_DRAWABLE, "drawable",  "Input drawable"   },
  162.     { GIMP_PDB_INT32,    "width",     "Width of output"  },
  163.     { GIMP_PDB_INT32,    "height",    "Height of output" }
  164.   };
  165.  
  166.   static GimpParamDef ret[] =
  167.   {
  168.     { GIMP_PDB_IMAGE,    "image_out", "Output image"     }
  169.   };
  170.  
  171.   /* use gimp_data_directory()/IBFNBQH-help for path */
  172.   help_path = g_build_filename (gimp_data_directory(), "IBFNBQH-help", NULL);
  173.   help_uri = g_filename_to_uri (help_path, NULL, NULL);
  174.   gimp_plugin_help_register (help_path, help_uri);
  175.   g_free (help_path);
  176.   g_free (help_uri);
  177.  
  178.   gimp_install_procedure (
  179.     "IBFNBQH",
  180.     "IBFNBQH",
  181. "Average matching image upsizing with Interpolatory Box Filtered Natural Biquadratic Histosplines.",
  182.     "Nicolas Robidoux, Adam Turcotte",
  183.     "Copyright Nicolas Robidoux, Adam Turcotte",
  184.     "2008",
  185.     "IBF_NBQH",
  186.     "RGB*, GRAY*",
  187.     GIMP_PLUGIN,
  188.     G_N_ELEMENTS (args),
  189.     G_N_ELEMENTS (ret),
  190.     args,
  191.     ret);
  192.  
  193.   gimp_plugin_menu_register ("IBFNBQH",
  194.                              "<Image>/Image");
  195.  
  196. }
  197.  
  198. static void
  199. run (const gchar      *name,
  200.      gint              nparams,
  201.      const GimpParam  *param,
  202.      gint             *nreturn_vals,
  203.      GimpParam       **return_vals)
  204. {
  205. #ifdef TIMING_REPORT
  206.   GTimer *timer;
  207. #endif
  208.  
  209.   static GimpParam  values[2];
  210.   GimpPDBStatusType status = GIMP_PDB_SUCCESS;
  211.   GimpRunMode       run_mode;
  212.   GimpDrawable     *drawable;
  213.   gint32            image_ID;
  214.   gint32            new_image_ID;
  215.  
  216.   gint m;        /* Height (in pixels) of the input image. */
  217.   gint n;        /* Width (in pixels) of the input image. */
  218.   gint mm;       /* Height (in pixels) of the output image. */
  219.   gint nn;       /* Width (in pixels) of the output image. */
  220.   gint channels; /* Number of colour channels of the input image. */
  221.  
  222.   /*
  223.    * Setting mandatory output values:
  224.    */
  225.   *nreturn_vals = 2;
  226.   *return_vals  = values;
  227.  
  228.   values[0].type = GIMP_PDB_STATUS;
  229.   values[0].data.d_status = status;
  230.  
  231.   /*
  232.    * Getting run_mode (we won't display a dialog if we are in
  233.    * NONINTERACTIVE mode):
  234.    */
  235.   run_mode = param[0].data.d_int32;
  236.  
  237.   image_ID = param[1].data.d_int32;
  238.  
  239.   /*  
  240.    * Get the specified drawable:  
  241.    */
  242.   drawable = gimp_drawable_get (param[2].data.d_drawable);
  243.  
  244.   _ibfnbqh_height = m = drawable->height;
  245.   _ibfnbqh_width = n = drawable->width;
  246.   channels = gimp_drawable_bpp (drawable->drawable_id);
  247.  
  248.   /*
  249.    * Display warning if image not at least 8x8
  250.    */
  251.   if ( MIN( n, m ) < 8 )
  252.   {
  253.     g_message ("Only enlarges images at least\n8 pixels wide and tall.\n");
  254.     return;
  255.   }
  256.  
  257.   switch (run_mode)
  258.   {
  259.     case GIMP_RUN_INTERACTIVE:
  260.       /* Initialize: */
  261.       svals.width  = n;
  262.       svals.height = m;
  263.       _ibfnbqh_m_over_n = (gdouble) m / n;
  264.  
  265.       /* Display the dialog: */
  266.       if (! scale_dialog (image_ID, m, n))
  267.         return;
  268.  
  269.       nn = svals.width;
  270.       mm = svals.height;
  271.  
  272.       /*
  273.        * Set up a tile cache large enough to contain a column of input
  274.        * tiles, or a row of output tiles.
  275.        */
  276.       gimp_tile_cache_ntiles (((nn>m) ? nn : m) / gimp_tile_width () + 1);
  277.  
  278. #ifdef TIMING_REPORT
  279.       timer = g_timer_new ();
  280. #ifdef EIGHT_BIT_STORAGE
  281.       g_print ("IBFNBQH (with uint8 coefficient storage) is taking...\n");
  282. #else
  283.       g_print ("IBFNBQH (with uint16 coefficient storage) is taking...\n");
  284. #endif
  285. #endif
  286.       new_image_ID = scale_up(image_ID, drawable, m, n, mm, nn, channels);
  287. #ifdef TIMING_REPORT
  288.       g_print ("%g seconds to complete the computation.\n",
  289.                g_timer_elapsed (timer, NULL));
  290.       g_timer_destroy (timer);
  291. #endif
  292.       values[1].type = GIMP_PDB_IMAGE;
  293.       values[1].data.d_image = new_image_ID;
  294.       gimp_displays_flush ();  /* flush to send data to core */
  295.       gimp_drawable_detach (drawable); /* then detach */
  296.  
  297.       /* Set options in the core */
  298.       gimp_set_data ("IBFNBQH", &svals, sizeof (ScaleVals));
  299.  
  300.       if (new_image_ID != -1)
  301.         gimp_display_new(new_image_ID);
  302.  
  303.       break;
  304.  
  305.     case GIMP_RUN_NONINTERACTIVE:
  306.       if (nparams != 5)
  307.         status = GIMP_PDB_CALLING_ERROR;
  308.       if (status == GIMP_PDB_SUCCESS)
  309.       {
  310.         svals.width = param[3].data.d_int32;
  311.         svals.height = param[4].data.d_int32;
  312.  
  313.         nn = svals.width;
  314.         mm = svals.height;
  315.  
  316.         /*
  317.          * Set up a tile cache large enough to contain a column of input
  318.          * tiles, or a row of output tiles.
  319.          */
  320.         gimp_tile_cache_ntiles (((nn>m) ? nn : m) / gimp_tile_width () + 1);
  321.  
  322. #ifdef TIMING_REPORT
  323.         timer = g_timer_new ();
  324. #ifdef EIGHT_BIT_STORAGE
  325.         g_print ("IBFNBQH (with int8 coefficient storage) is taking...\n");
  326. #else
  327.         g_print ("IBFNBQH (with int16 coefficient storage) is taking...\n");
  328. #endif
  329. #endif
  330.         new_image_ID = scale_up(image_ID, drawable, m, n, mm, nn, channels);
  331. #ifdef TIMING_REPORT
  332.         g_print ("%g seconds to complete the computation.\n",
  333.                  g_timer_elapsed (timer, NULL));
  334.         g_timer_destroy (timer);
  335. #endif
  336.         values[1].type = GIMP_PDB_IMAGE;
  337.         values[1].data.d_image = new_image_ID;
  338.         gimp_displays_flush ();  /* flush to send data to core */
  339.         gimp_drawable_detach (drawable); /* then detach */
  340.  
  341.         /* Set options in the core */
  342.         gimp_set_data ("IBFNBQH", &svals, sizeof (ScaleVals));
  343.       }
  344.       break;
  345.  
  346.     case GIMP_RUN_WITH_LAST_VALS:
  347.       gimp_get_data ("IBFNBQH", &svals);
  348.       nn = svals.width;
  349.       mm = svals.height;
  350.  
  351.       /*
  352.        * Set up a tile cache large enough to contain a column of input
  353.        * tiles, or a row of output tiles.
  354.        */
  355.       gimp_tile_cache_ntiles (((nn>m) ? nn : m) / gimp_tile_width () + 1);
  356.  
  357. #ifdef TIMING_REPORT
  358.       timer = g_timer_new ();
  359. #ifdef EIGHT_BIT_STORAGE
  360.       g_print ("IBFNBQH (with uint8 coefficient storage) is taking...\n");
  361. #else
  362.       g_print ("IBFNBQH (with uint16 coefficient storage) is taking...\n");
  363. #endif
  364. #endif
  365.       new_image_ID = scale_up(image_ID, drawable, m, n, mm, nn, channels);
  366. #ifdef TIMING_REPORT
  367.       g_print ("%g seconds to complete the computation.\n",
  368.                g_timer_elapsed (timer, NULL));
  369.       g_timer_destroy (timer);
  370. #endif
  371.       values[1].type = GIMP_PDB_IMAGE;
  372.       values[1].data.d_image = new_image_ID;
  373.       gimp_displays_flush ();  /* flush to send data to core */
  374.       gimp_drawable_detach (drawable); /* then detach */
  375.  
  376.       break;
  377.  
  378.     default:
  379.       break;
  380.   }
  381.  
  382.   return;
  383. }
  384.  
  385. static gint32
  386. scale_up (gint32        image_ID,
  387.           GimpDrawable *drawable,
  388.           gint          m,
  389.           gint          n,
  390.           gint          mm,
  391.           gint          nn,
  392.           gint          channels)
  393. {
  394.   gint first_past_ii[m-1];
  395.   gint first_past_jj[n-1];
  396.  
  397.   gfloat top[mm];
  398.   gfloat middle[mm];
  399.   gfloat bottom[mm];
  400.   gfloat farbottom[mm];
  401.  
  402.   gfloat left[nn];
  403.   gfloat center[nn];
  404.   gfloat right[nn];
  405.   gfloat farright[nn];
  406.  
  407.   GimpDrawable *new_drawable;
  408.   gint32 new_image_ID;
  409.   gint32 new_layer_ID;
  410.  
  411.   GimpPixelRgn input_image, output_image;
  412.  
  413.   /*
  414.    * Useful iterator bounds:
  415.    */
  416.   gint nn_times_channels                = nn * channels;
  417.   gint n_times_channels                 = n * channels;
  418.   gint n_minus_1_times_channels         = (n-1) * channels;
  419.   gint n_minus_7_times_channels         = (n-7) * channels;
  420.   gint n_minus_8_times_channels         = n_minus_7_times_channels - channels;
  421.   gint n_minus_2                        = n - 2;
  422.   gint m_times_channels                 = m * channels;
  423.   gint m_minus_6_times_channels         = (m - 6) * channels;
  424.   gint m_minus_7_times_channels         = m_minus_6_times_channels - channels;
  425.   gint m_minus_2                        = m - 2;
  426.  
  427.   /*
  428.    * Utility iterators:
  429.    */
  430.   gint k; /* Index related to the number of pixel values in an image
  431.              row (or column). If the image has, say, three channels,
  432.              to traverse a row one needs to visit 3n entries; k is
  433.              used for that. */
  434.   gint c; /* Index (usually) related to the number of color channels
  435.              (3 for RGB, 1 for Grayscale...). */
  436.  
  437.   gint kk;
  438.  
  439.   /*
  440.    * Row and column indices:
  441.    */
  442.   gint i; /* Index of the input row under consideration. */
  443.   gint j; /* Index of the input column under consideration. */
  444.  
  445.   gint ii; /* Index of the output row under consideration. */
  446.   gint last_ii; /* Local loop bound for ii. */
  447.  
  448.   gint jj; /* Index of the output column under consideration. */
  449.   gint last_jj; /* Local loop bound for jj. */
  450.  
  451.   /*
  452.    * Computed coefficients/partial values:
  453.    */
  454.   gfloat top_ii;
  455.   gfloat middle_ii;
  456.   gfloat bottom_ii;
  457.   gfloat farbottom_ii;
  458.  
  459.   gfloat left_jj;
  460.   gfloat center_jj;
  461.   gfloat right_jj;
  462.   gfloat farright_jj;
  463.  
  464. #ifdef TIMING_REPORT
  465.   GTimer *local_timer;
  466.   gfloat local_time;
  467. #endif
  468.  
  469.   gfloat a_top_left[channels];
  470.   gfloat a_top_center[channels];
  471.   gfloat a_top_right[channels];
  472.   gfloat a_top_farright[channels];
  473.  
  474.   gfloat a_middle_left[channels];
  475.   gfloat a_middle_center[channels];
  476.   gfloat a_middle_right[channels];
  477.   gfloat a_middle_farright[channels];
  478.  
  479.   gfloat a_bottom_left[channels];
  480.   gfloat a_bottom_center[channels];
  481.   gfloat a_bottom_right[channels];
  482.   gfloat a_bottom_farright[channels];
  483.  
  484.   gfloat a_farbottom_left[channels];
  485.   gfloat a_farbottom_center[channels];
  486.   gfloat a_farbottom_right[channels];
  487.   gfloat a_farbottom_farright[channels];
  488.  
  489.   gfloat coef_left[channels];
  490.   gfloat coef_center[channels];
  491.   gfloat coef_right[channels];
  492.   gfloat coef_farright[channels];
  493.  
  494.   /*
  495.    * Computation of the B-spline coefficients:
  496.    *
  497.    * For natural boundary conditions:
  498.    * C0      = 1/5                
  499.    * C1      = 5/19                = 1 / ( 4 - C0 )
  500.    * C2      = 19/71               = 1 / ( 4 - C1 )
  501.    * C3      = 71/265              = 1 / ( 4 - C2 )
  502.    * C4      = 265/989             = 1 / ( 4 - C3 )
  503.    * C5      = 989/3691            = 1 / ( 4 - C4 )
  504.    * CINFTY  = 2 - sqrt(3)         = limit of this recurrence relation
  505.    * CLAST   = 1 / ( 5 - C_INFTY ) = ( 3 - sqrt(3) ) / 6
  506.    *
  507.    * Convergence of the coefficients for natural boundary conditions:
  508.    * In 32 bit single precision, C6 and all subsequent ones are equal
  509.    * to CINFTY. Although C5 is different from CINFTY in this
  510.    * precision, it differs from it only by 1.271371644 * 10^-7, just
  511.    * enough to be different in standard single precision, but
  512.    * nonetheless an insignificant amount (.2679491f instead of
  513.    * .2679492f). However, in the interest of preventing integer
  514.    * overflow, we actually use C5 as a proxy for CINFTY in the
  515.    * computations which take place prior to storing partially computed
  516.    * results into int8s. This is discussed further below.
  517.    */
  518. #define C0_F        .2000000f
  519. #define MINUS_C0_F -.2000000f
  520. #define C1_F        .2631579f
  521. #define MINUS_C1_F -.2631579f
  522. #define C2_F        .2676056f
  523. #define MINUS_C2_F -.2676056f
  524. #define C3_F        .2679245f
  525. #define MINUS_C3_F -.2679245f
  526. #define C4_F        .2679474f
  527. #define MINUS_C4_F -.2679474f
  528. #define C5_F        .2679491f
  529. #define MINUS_C5_F -.2679491f
  530. #define CLAST_F     .2113249f
  531.   /*
  532.    * So as to sure that we don't overflow when converting partially
  533.    * computed coefficients to integer data types (without explicitly
  534.    * clamping), in the computations which come prior to the cast we
  535.    * use C5 as a proxy for CINFTY: C5 turns out to be the truncated
  536.    * value of CINFTY, which is what we want. In the computations which
  537.    * take place after the partially computed B-spline coefficients
  538.    * have been converted back to floats, we use the correctly rounded
  539.    * value of CINFTY.
  540.    */
  541. #define CINFTY_F        .2679492f
  542. #define MINUS_CINFTY_F -.2679492f
  543.  
  544.   gfloat *a_row_p1, *a_row_p2;
  545.  
  546.   /*
  547.    * Rescaled partially computed B-Spline coefficients are stored in
  548.    * a.
  549.    */
  550. #ifdef EIGHT_BIT_STORAGE
  551.   guchar *a, *a_ptr;
  552.   a = g_new (guchar,  n * m_times_channels);
  553. #else
  554.   guint16 *a, *a_ptr;
  555.   a = g_new (guint16, n * m_times_channels);
  556. #endif
  557.  
  558. #ifdef TIMING_REPORT
  559.   local_timer = g_timer_new ();
  560.   local_time  = g_timer_elapsed (local_timer, NULL);
  561. #endif
  562.  
  563.   first_past_index(m, mm, first_past_ii);
  564.   first_past_index(n, nn, first_past_jj);
  565.  
  566. #ifdef TIMING_REPORT
  567.   g_print ("%g seconds for the first_past_index calls.\n",
  568.            g_timer_elapsed (local_timer, NULL)-local_time);
  569.   local_time = g_timer_elapsed (local_timer, NULL);
  570. #endif
  571.  
  572.   /*
  573.    * Computation of the tensor components of the linear transformation
  574.    * from B-spline coefficents to fine cell averages (actually fine
  575.    * cell integrals, since we have alread rescaled by the reciprocal
  576.    * of the fine cell areas):
  577.    */
  578.   smooth_coarse_to_fine_coefficients(m, mm, first_past_ii,
  579.                                      top, middle, bottom, farbottom);
  580.   smooth_coarse_to_fine_coefficients(n, nn, first_past_jj,
  581.                                      left, center, right, farright);
  582.  
  583. #ifdef TIMING_REPORT
  584.   g_print ("%g seconds for the smooth_coarse_to_fine_coefficients calls.\n",
  585.            g_timer_elapsed (local_timer, NULL)-local_time);
  586. #endif
  587.  
  588.   /*
  589.    * Initialize the input pixel region:
  590.    */
  591.   gimp_pixel_rgn_init (&input_image, drawable, 0, 0, n, m, FALSE, FALSE);
  592.  
  593.   /*
  594.    * Set up an informative progress bar:
  595.    */
  596.   gimp_progress_init (
  597.     "Upsizing with average matching biquadratic histosplines");
  598.  
  599.   gimp_progress_update(.005);
  600.  
  601. #ifdef TIMING_REPORT
  602.   local_time  = g_timer_elapsed (local_timer, NULL);
  603. #endif
  604.  
  605.   {
  606.     /*
  607.      * Column-based (row by row) forward substitution, performed one
  608.      * column at a time.
  609.      */    
  610.  
  611.     guchar *input_col, *input_col_p; /* Storage and moving pointer for
  612.                                         input image columns. */
  613.     gfloat *a_col, *a_col_p1, *a_col_p2; /* The computation is
  614.                                             performed in the a_col
  615.                                             array, accessed using the
  616.                                             two corresponding
  617.                                             pointers. */
  618.  
  619.     input_col = g_new (guchar, m_times_channels);    
  620.     a_col     = g_new (gfloat, m_times_channels);
  621.  
  622.     for (j=0; j<n; j++)
  623.     {
  624.       /*
  625.        * Implicitly, the code rescales the 0..255 uchar values so that
  626.        * they fit in -255 to 255. Because the infinity norm of the
  627.        * inverse of the tensor component of the "interpolation matrix"
  628.        * is 1/2, the computed coefficients fit into -127.5 to 127.5,
  629.        * which we later encode into 0 to 65535.
  630.        *
  631.        * The purpose of FIT_TO_RANGE is to turn pixel values ranging
  632.        * from 0 to 255 to values ranging from -255 to 255. Values in
  633.        * this range are mapped (arbitrarily close to when the matrix
  634.        * is large) values strictly in the range -127.5 to 127.5 by the
  635.        * inverse of the matrix A =
  636.        *
  637.        * 5 1
  638.        * 1 4 1
  639.        *   1 4 1
  640.        *     1 4 1
  641.        *       1 4 1
  642.        *         ...
  643.        *           1 5
  644.        *
  645.        * which is the matrix for natural boundary conditions.
  646.        *
  647.        * FIT_TO_SCALE scales by a factor of 2 (and shifts).
  648.        *
  649.        * One can understand the factor of two as being part of solving
  650.        * when the matrix is actually half the one shown above. This
  651.        * half-matrix has an interpretation in terms of rescaling
  652.        * B-splines.
  653.        *
  654.        * The advantage of this interpretation is that this scaling
  655.        * (asymptotically) keeps the range invariant:
  656.        *
  657.        * If all entries of y are within -127.5 to 127.5 (inclusive),
  658.        * and x is the solution of A'x = y, then all entries of x are
  659.        * within the same range (exclusive), provided
  660.        *
  661.        * A' =
  662.        * 2.5 .5
  663.        * .5   2 .5
  664.        *     .5  2 .5
  665.        *        .5  2 .5
  666.        *            ...
  667.        *           .5 2.5
  668.        *
  669.        * Halving the matrix is equivalent to doubling the RHS.
  670.        *
  671.        * Yet a third interpretation: In order to store results
  672.        * using the first version of the matrix, we double them (at
  673.        * any stage), and then have to halve them when we want to
  674.        * use them.
  675.        *
  676.        * The important thing to keep in mind for later is that the
  677.        * values are scaled by 2.
  678.        */
  679.  
  680.       /*
  681.        * ADAM AND NICOLAS: When we have time, this should be
  682.        * implemented as a (implicit) cast, a bit shift, and an integer
  683.        * add (then, no flop is involved here):
  684.        */
  685. #define FIT_TO_RANGE(x) ( 2 * ( (int) (x) ) - 255 )
  686.       /*
  687.        * The above can be rewritten using integer arithmetic, at the
  688.        * cost of an additional (explicit) intermediate cast, as
  689.        * follows. This does not seem to speed up the code much.
  690.        */
  691.       /* #define FIT_TO_RANGE(x) ({gint in = (x); (in+in) - 255;}) */
  692.  
  693.       /*
  694.        * The code for this section has been structured so as to make
  695.        * obvious multiply-add operations. Whether this is exploited by
  696.        * the compiler is another story. Currently, gcc is not C99
  697.        * compliant: for this reason, the current benefit of this is
  698.        * probably nil. It could, however, make a significant
  699.        * difference with another compiler.
  700.        */
  701.  
  702.       /*
  703.        * Read one column of intensity values into input_col.
  704.        */
  705.       gimp_pixel_rgn_get_col (&input_image, input_col, j, 0, m);
  706.  
  707.       /*
  708.        * Non-asymptotic LU factorization entries:
  709.        */
  710.       input_col_p = input_col;
  711.       a_col_p1    = a_col;
  712.       a_col_p2    = a_col;
  713.      
  714.       for (c=channels; c; c--, a_col_p1++, input_col_p++)
  715.         *a_col_p1 =                        FIT_TO_RANGE( *input_col_p );
  716.  
  717.       for (c=channels; c; c--, a_col_p1++, a_col_p2++, input_col_p++)
  718.         *a_col_p1 = *a_col_p2 * MINUS_C0_F + FIT_TO_RANGE( *input_col_p );
  719.  
  720.       for (c=channels; c; c--, a_col_p1++, a_col_p2++, input_col_p++)
  721.         *a_col_p1 = *a_col_p2 * MINUS_C1_F + FIT_TO_RANGE( *input_col_p );
  722.  
  723.       for (c=channels; c; c--, a_col_p1++, a_col_p2++, input_col_p++)
  724.         *a_col_p1 = *a_col_p2 * MINUS_C2_F + FIT_TO_RANGE( *input_col_p );
  725.  
  726.       for (c=channels; c; c--, a_col_p1++, a_col_p2++, input_col_p++)
  727.         *a_col_p1 = *a_col_p2 * MINUS_C3_F + FIT_TO_RANGE( *input_col_p );
  728.  
  729.       for (c=channels; c; c--, a_col_p1++, a_col_p2++, input_col_p++)
  730.         *a_col_p1 = *a_col_p2 * MINUS_C4_F + FIT_TO_RANGE( *input_col_p );
  731.  
  732.       /*
  733.        * This is the first spot where we use C5 as a proxy for CINFTY.
  734.        */
  735.       /*
  736.        * Asymptotic (within roundoff) LU factorization entries, except
  737.        * for the very last one:
  738.        */
  739.       for (k=m_minus_7_times_channels; k; k--,
  740.              a_col_p1++, a_col_p2++, input_col_p++)
  741.         *a_col_p1 = *a_col_p2 * MINUS_C5_F + FIT_TO_RANGE( *input_col_p );
  742.  
  743.       /*
  744.        * Last step of column-based (row by row) forward substitution,
  745.        * and first step of column-based backward substitution,
  746.        * performed in one step on the very last entries:
  747.        */
  748.       for (c=channels; c; c--, a_col_p1++, a_col_p2++, input_col_p++)
  749.         *a_col_p1 = ( *a_col_p2 * MINUS_CINFTY_F
  750.                       + FIT_TO_RANGE( *input_col_p ) ) * CLAST_F;
  751.  
  752.       /*
  753.        * Remainder of the row-based back substitution (we just took
  754.        * care of the very last row, which is why the indices are
  755.        * rewound the way they are).
  756.        */
  757.       a_col_p1--;
  758.       a_col_p2 = a_col_p1;
  759.       for (c=channels; c; c--)
  760.         a_col_p1--;
  761.  
  762.       /*
  763.        * This is the second spot where we use C5 as a proxy for
  764.        * CINFTY:
  765.        */
  766.       for (k=m_minus_6_times_channels; k; k--, a_col_p1--, a_col_p2--)
  767.         *a_col_p1 = ( *a_col_p1 - *a_col_p2 ) * C5_F;
  768.  
  769.       for (c=channels; c; c--, a_col_p1--, a_col_p2--)
  770.         *a_col_p1 = ( *a_col_p1 - *a_col_p2 ) * C4_F;
  771.  
  772.       for (c=channels; c; c--, a_col_p1--, a_col_p2--)
  773.         *a_col_p1 = ( *a_col_p1 - *a_col_p2 ) * C3_F;
  774.  
  775.       for (c=channels; c; c--, a_col_p1--, a_col_p2--)
  776.         *a_col_p1 = ( *a_col_p1 - *a_col_p2 ) * C2_F;
  777.  
  778.       for (c=channels; c; c--, a_col_p1--, a_col_p2--)
  779.         *a_col_p1 = ( *a_col_p1 - *a_col_p2 ) * C1_F;
  780.  
  781.       for (c=channels; c; c--, a_col_p1--, a_col_p2--)
  782.         *a_col_p1 = ( *a_col_p1 - *a_col_p2 ) * C0_F;
  783.  
  784.       a_col_p1 = a_col;
  785.       a_ptr = a+j*channels;
  786.  
  787.       /*
  788.        * Pack the partially computed coefficients into uint16/uint8:
  789.        *
  790.        * Convert the column of partially computed B-spline
  791.        * coefficients into uint8 or uint16, and append the column to
  792.        * the end of uint array of B-spline coefficients a.
  793.        *
  794.        * In the case of uint16, we store the range -127.5 to 127.5 by
  795.        * affine transforming it to 0 to 65535, combined with rounding
  796.        * (toward +infty for no reason other than expediency; Banker's
  797.        * rounding would be fine). Using cast (= truncate) to effect
  798.        * the rounding, we get that the rounded value is obtained by
  799.        * multiplying by 257 to get values in the range -32767.5 to
  800.        * 32676.5, then adding 32768 so that truncation rounds
  801.        * correctly (with ties resolved towards +infinity
  802.        * ("white"---definitely not Banker's rounding)). Note that the
  803.        * inverse of A keeps values strictly within -127.5 to 127.5
  804.        * (these bounds are NOT included).
  805.        *
  806.        * In the case of storage into uint8, things are a bit
  807.        * simpler. However, we must still perform the shift in float
  808.        * arithmetic prior to the implicit cast, because otherwise we
  809.        * mess up the emulation of rounding. Rounding emulation is the
  810.        * reason we use uint8 even though int8 are naturally in a range
  811.        * which most closely resembles -127.5 to 127.5.
  812.        */
  813.       for (i=m; i; i--, a_ptr+=n_minus_1_times_channels)
  814.         for (c=channels; c; c--, a_col_p1++, a_ptr++)
  815. #ifdef EIGHT_BIT_STORAGE
  816.           *a_ptr = *a_col_p1 + 128.f;
  817. #else
  818.           *a_ptr = 257.f * *a_col_p1 + 32768.f;
  819. #endif
  820.     }
  821.  
  822. #ifdef TIMING_REPORT
  823.   g_print ("%g seconds for forward substitution.\n",
  824.            g_timer_elapsed (local_timer, NULL)-local_time);
  825.   g_timer_destroy (local_timer);
  826. #endif
  827.  
  828.     gimp_progress_update(.005+((gdouble)m*n)/((gdouble)mm*nn)*.14);
  829.  
  830.     g_free(a_col);
  831.     g_free(input_col);
  832.   }
  833.  
  834.   /*
  835.    * Create new image to place scaled drawable inside.
  836.    */
  837.   new_image_ID = gimp_image_new (nn, mm, gimp_image_base_type (image_ID));
  838.   new_layer_ID = gimp_layer_new (new_image_ID, "Background",
  839.                                  nn, mm,
  840.                                  gimp_drawable_type (drawable->drawable_id),
  841.                                  100,
  842.                                  GIMP_NORMAL_MODE);
  843.   gimp_image_add_layer (new_image_ID, new_layer_ID, 0);
  844.   new_drawable = gimp_drawable_get (new_layer_ID);
  845.  
  846.   /*
  847.    * Init the output pixel region:
  848.    */
  849.   gimp_pixel_rgn_init (&output_image, new_drawable, 0, 0, nn, mm, TRUE, TRUE);
  850.  
  851.   {
  852.     gfloat *a_middle, *a_bottom, *a_farbottom, *a_top;
  853.     gfloat *a_middle_p, *a_bottom_p, *a_farbottom_p, *a_top_p, *a_temp;
  854.    
  855.     /*
  856.      * Storage for a pixel row to be output to new image buffer:
  857.      */
  858.     guchar *output_row;
  859.     guchar *output_row_ptr;
  860.  
  861.     gint nn_minus_1 = nn-1;
  862.  
  863.     a_top       = g_new (gfloat, n_times_channels);
  864.     a_middle    = g_new (gfloat, n_times_channels);
  865.     a_bottom    = g_new (gfloat, n_times_channels);
  866.     a_farbottom = g_new (gfloat, n_times_channels);
  867.  
  868.     /*
  869.      * Temporary storage for a row of the output image:
  870.      */
  871.     output_row  = g_new (guchar, nn_times_channels);
  872.  
  873.     a_ptr = a;
  874.  
  875.     /*
  876.      * ADAM AND NICOLAS: WHEN WE HAVE TIME, REWRITE SO NO FLOP IS
  877.      * INVOLVED (using integer arithmetic):
  878.      */
  879. #ifdef EIGHT_BIT_STORAGE
  880. #define UNSCALING ( (float) (.5) )
  881. #define UNSCALED_SHIFT (-85)
  882. #else
  883. #define UNSCALING ( (float) (1./514.) )
  884. #define UNSCALED_SHIFT (-21845)
  885. #endif
  886.  
  887.     /*
  888.      * Row solve on a_middle:
  889.      */
  890.  
  891.     a_row_p1 = a_middle;
  892.  
  893.     for (kk=n_times_channels; kk; kk--, a_row_p1++, a_ptr++)
  894.       *a_row_p1 = ( (gint) *a_ptr + UNSCALED_SHIFT ) * UNSCALING;
  895.  
  896.     /*
  897.      * Forward substitution:
  898.      */
  899.  
  900.     /*
  901.      * We'll take care of the very first pixel on the way back.
  902.      */
  903.     a_row_p1 = a_middle+channels;
  904.     a_row_p2 = a_middle;
  905.  
  906.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  907.       *a_row_p1 += *a_row_p2 * MINUS_C0_F;
  908.  
  909.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  910.       *a_row_p1 += *a_row_p2 * MINUS_C1_F;
  911.  
  912.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  913.       *a_row_p1 += *a_row_p2 * MINUS_C2_F;
  914.  
  915.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  916.       *a_row_p1 += *a_row_p2 * MINUS_C3_F;
  917.  
  918.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  919.       *a_row_p1 += *a_row_p2 * MINUS_C4_F;
  920.  
  921.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  922.       *a_row_p1 += *a_row_p2 * MINUS_C5_F;
  923.  
  924.     /*
  925.      * Because integer overflow is not a risk any more---we do clamp
  926.      * values later on, and we'd have to even in exact arithmetic---we
  927.      * use the rounded value of CINFTY instead of its truncated one.
  928.      */
  929.     for (kk=n_minus_8_times_channels; kk; kk--, a_row_p1++, a_row_p2++)
  930.       *a_row_p1 += *a_row_p2 * MINUS_CINFTY_F;
  931.  
  932.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  933.       *a_row_p1 = ( *a_row_p2 * MINUS_CINFTY_F + *a_row_p1 ) * CLAST_F;
  934.  
  935.     a_row_p1--;
  936.     a_row_p2 = a_row_p1;
  937.     for (c=channels; c; c--)
  938.       a_row_p1--;
  939.  
  940.     for (kk=n_minus_7_times_channels; kk; kk--, a_row_p1--, a_row_p2--)
  941.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * CINFTY_F;
  942.  
  943.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  944.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C5_F;
  945.  
  946.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  947.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C4_F;
  948.  
  949.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  950.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C3_F;
  951.  
  952.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  953.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C2_F;
  954.  
  955.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  956.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C1_F;
  957.  
  958.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  959.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C0_F;
  960.  
  961.     /*
  962.      * Row solve on a_bottom:
  963.      */
  964.  
  965.     a_row_p1 = a_bottom;
  966.  
  967.     for (kk=n_times_channels; kk; kk--, a_row_p1++, a_ptr++)
  968.       *a_row_p1 = ( (gint) *a_ptr + UNSCALED_SHIFT ) * UNSCALING;
  969.  
  970.     /*
  971.      * Forward substitution:
  972.      */
  973.  
  974.     /*
  975.      * We'll take care of the very first pixel on the way back.
  976.      */
  977.     a_row_p1 = a_bottom+channels;
  978.     a_row_p2 = a_bottom;
  979.  
  980.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  981.       *a_row_p1 += *a_row_p2 * MINUS_C0_F;
  982.  
  983.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  984.       *a_row_p1 += *a_row_p2 * MINUS_C1_F;
  985.  
  986.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  987.       *a_row_p1 += *a_row_p2 * MINUS_C2_F;
  988.  
  989.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  990.       *a_row_p1 += *a_row_p2 * MINUS_C3_F;
  991.  
  992.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  993.       *a_row_p1 += *a_row_p2 * MINUS_C4_F;
  994.  
  995.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  996.       *a_row_p1 += *a_row_p2 * MINUS_C5_F;
  997.  
  998.     /*
  999.      * Because integer overflow is not a risk any more---we do clamp
  1000.      * values later on, and we'd have to even in exact arithmetic---we
  1001.      * use the rounded value of CINFTY instead of its truncated one.
  1002.      */
  1003.     for (kk=n_minus_8_times_channels; kk; kk--, a_row_p1++, a_row_p2++)
  1004.       *a_row_p1 += *a_row_p2 * MINUS_CINFTY_F;
  1005.  
  1006.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1007.       *a_row_p1 = ( *a_row_p2 * MINUS_CINFTY_F + *a_row_p1 ) * CLAST_F;
  1008.  
  1009.     a_row_p1--;
  1010.     a_row_p2 = a_row_p1;
  1011.     for (c=channels; c; c--)
  1012.       a_row_p1--;
  1013.  
  1014.     for (kk=n_minus_7_times_channels; kk; kk--, a_row_p1--, a_row_p2--)
  1015.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * CINFTY_F;
  1016.  
  1017.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1018.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C5_F;
  1019.  
  1020.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1021.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C4_F;
  1022.  
  1023.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1024.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C3_F;
  1025.  
  1026.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1027.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C2_F;
  1028.  
  1029.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1030.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C1_F;
  1031.  
  1032.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1033.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C0_F;
  1034.  
  1035.     /*
  1036.      * Row solve on a_farbottom:
  1037.      */
  1038.  
  1039.     a_row_p1 = a_farbottom;
  1040.  
  1041.     for (kk=n_times_channels; kk; kk--, a_row_p1++, a_ptr++)
  1042.       *a_row_p1 = ( (gint) *a_ptr + UNSCALED_SHIFT ) * UNSCALING;
  1043.  
  1044.     /*
  1045.      * Forward substitution:
  1046.      */
  1047.  
  1048.     /*
  1049.      * We'll take care of the very first pixel on the way back.
  1050.      */
  1051.     a_row_p1 = a_farbottom+channels;
  1052.     a_row_p2 = a_farbottom;
  1053.  
  1054.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1055.       *a_row_p1 += *a_row_p2 * MINUS_C0_F;
  1056.  
  1057.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1058.       *a_row_p1 += *a_row_p2 * MINUS_C1_F;
  1059.  
  1060.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1061.       *a_row_p1 += *a_row_p2 * MINUS_C2_F;
  1062.  
  1063.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1064.       *a_row_p1 += *a_row_p2 * MINUS_C3_F;
  1065.  
  1066.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1067.       *a_row_p1 += *a_row_p2 * MINUS_C4_F;
  1068.  
  1069.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1070.       *a_row_p1 += *a_row_p2 * MINUS_C5_F;
  1071.  
  1072.     /*
  1073.      * Because integer overflow is not a risk any more---we do clamp
  1074.      * values later on, and we'd have to even in exact arithmetic---we
  1075.      * use the rounded value of CINFTY instead of its truncated one.
  1076.      */
  1077.     for (kk=n_minus_8_times_channels; kk; kk--, a_row_p1++, a_row_p2++)
  1078.       *a_row_p1 += *a_row_p2 * MINUS_CINFTY_F;
  1079.  
  1080.     for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1081.       *a_row_p1 = ( *a_row_p2 * MINUS_CINFTY_F + *a_row_p1 ) * CLAST_F;
  1082.  
  1083.     a_row_p1--;
  1084.     a_row_p2 = a_row_p1;
  1085.     for (c=channels; c; c--)
  1086.       a_row_p1--;
  1087.  
  1088.     for (kk=n_minus_7_times_channels; kk; kk--, a_row_p1--, a_row_p2--)
  1089.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * CINFTY_F;
  1090.  
  1091.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1092.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C5_F;
  1093.  
  1094.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1095.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C4_F;
  1096.  
  1097.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1098.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C3_F;
  1099.  
  1100.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1101.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C2_F;
  1102.  
  1103.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1104.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C1_F;
  1105.  
  1106.     for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1107.       *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C0_F;
  1108.  
  1109.     gimp_progress_update(.005+((gdouble)m*n)/((gdouble)mm*nn)*.15);
  1110.  
  1111.     last_ii = first_past_ii[0];
  1112.    
  1113.     for (ii=0; ii<last_ii; ii++)
  1114.       {
  1115.         /*
  1116.          * Move output_row_ptr to the beginning of output_row:
  1117.          */
  1118.         output_row_ptr = output_row;
  1119.  
  1120.         middle_ii    = middle[ii];
  1121.         bottom_ii    = bottom[ii];
  1122.         farbottom_ii = farbottom[ii];
  1123.  
  1124.         a_middle_p    = a_middle;
  1125.         a_bottom_p    = a_bottom;
  1126.         a_farbottom_p = a_farbottom;
  1127.  
  1128.         for (c=0; c<channels; c++)
  1129.           {
  1130.             a_middle_center[c]      = *a_middle_p++;
  1131.             a_bottom_center[c]      = *a_bottom_p++;
  1132.             a_farbottom_center[c]   = *a_farbottom_p++;
  1133.           }
  1134.  
  1135.         for (c=0; c<channels; c++)
  1136.           {
  1137.             a_middle_right[c]       = *a_middle_p++;
  1138.             a_bottom_right[c]       = *a_bottom_p++;
  1139.             a_farbottom_right[c]    = *a_farbottom_p++;
  1140.           }
  1141.        
  1142.         for (c=0; c<channels; c++)
  1143.           {
  1144.             a_middle_farright[c]    = *a_middle_p++;
  1145.             a_bottom_farright[c]    = *a_bottom_p++;
  1146.             a_farbottom_farright[c] = *a_farbottom_p++;
  1147.           }        
  1148.  
  1149.         for (c=0; c<channels; c++)
  1150.           {
  1151.             coef_center[c]   =   middle_ii    *      a_middle_center[c]
  1152.                                + bottom_ii    *      a_bottom_center[c]
  1153.                                + farbottom_ii *   a_farbottom_center[c];
  1154.             coef_right[c]    =   middle_ii    *       a_middle_right[c]
  1155.                                + bottom_ii    *       a_bottom_right[c]
  1156.                                + farbottom_ii *    a_farbottom_right[c];
  1157.             coef_farright[c] =   middle_ii    *    a_middle_farright[c]
  1158.                                + bottom_ii    *    a_bottom_farright[c]
  1159.                                + farbottom_ii * a_farbottom_farright[c];
  1160.           }
  1161.  
  1162.         last_jj = first_past_jj[0];
  1163.        
  1164.         for (jj=0; jj<last_jj; jj++)
  1165.           {
  1166.             center_jj   =   center[jj];
  1167.             right_jj    =    right[jj];
  1168.             farright_jj = farright[jj];
  1169.  
  1170.             /*
  1171.              * F_ROUND_AND_CLAMP_0255(x) rounds x to the nearest integer,
  1172.              * clamping to values between 0 and 255, and casting it to a
  1173.              * guchar.
  1174.              *
  1175.              * The ".5f +" is so that truncation (which the cast from
  1176.              * float to int effects for positive values) performs
  1177.              * round to nearest. Negative values are clamped later on
  1178.              * to 0, so the fact that truncation brings things "up"
  1179.              * for negative results is not a concern. Halfway points
  1180.              * between nonnnegative integers are mapped to the larger
  1181.              * integer, which is acceptable---maybe even good if the
  1182.              * image is overexposed---but not important: round to
  1183.              * nearest, with any convention (banker's, toward
  1184.              * zero---which may have advantages---or toward +infty,
  1185.              * which is the case now, would also be fine).
  1186.              */
  1187. #define F_ROUND_AND_CLAMP_0255(x) ({ gint out = .5f + (x); (out > 255) ? (guchar)255 : ((out < 0) ? (guchar)0 : (guchar)out); })
  1188.  
  1189.             for (c=0; c<channels; c++)
  1190.               *output_row_ptr++ =
  1191.                 F_ROUND_AND_CLAMP_0255(   coef_center[c]   * center_jj
  1192.                                         + coef_right[c]    * right_jj
  1193.                                         + coef_farright[c] * farright_jj);
  1194.           }
  1195.  
  1196.         for (j=1; j<n_minus_2; j++)
  1197.           {
  1198.             for (c=0; c<channels; c++)
  1199.               {
  1200.                 a_middle_left[c]        =  a_middle_center[c];
  1201.                 a_middle_center[c]      =  a_middle_right[c];
  1202.                 a_middle_right[c]       =  a_middle_farright[c];
  1203.                 a_middle_farright[c]    = *a_middle_p++;
  1204.  
  1205.                 a_bottom_left[c]        =  a_bottom_center[c];
  1206.                 a_bottom_center[c]      =  a_bottom_right[c];
  1207.                 a_bottom_right[c]       =  a_bottom_farright[c];
  1208.                 a_bottom_farright[c]    = *a_bottom_p++;
  1209.  
  1210.                 a_farbottom_left[c]     =  a_farbottom_center[c];
  1211.                 a_farbottom_center[c]   =  a_farbottom_right[c];
  1212.                 a_farbottom_right[c]    =  a_farbottom_farright[c];
  1213.                 a_farbottom_farright[c] = *a_farbottom_p++;
  1214.  
  1215.                 coef_left[c]     =   middle_ii    *        a_middle_left[c]
  1216.                                    + bottom_ii    *        a_bottom_left[c]
  1217.                                    + farbottom_ii *     a_farbottom_left[c];
  1218.                 coef_center[c]   =   middle_ii    *      a_middle_center[c]
  1219.                                    + bottom_ii    *      a_bottom_center[c]
  1220.                                    + farbottom_ii *   a_farbottom_center[c];
  1221.                 coef_right[c]    =   middle_ii    *       a_middle_right[c]
  1222.                                    + bottom_ii    *       a_bottom_right[c]
  1223.                                    + farbottom_ii *    a_farbottom_right[c];
  1224.                 coef_farright[c] =   middle_ii    *    a_middle_farright[c]
  1225.                                    + bottom_ii    *    a_bottom_farright[c]
  1226.                                    + farbottom_ii * a_farbottom_farright[c];
  1227.               }
  1228.  
  1229.             last_jj = first_past_jj[j];
  1230.  
  1231.             for (; jj<last_jj; jj++)
  1232.               {
  1233.                 left_jj     =     left[jj];
  1234.                 center_jj   =   center[jj];
  1235.                 right_jj    =    right[jj];
  1236.                 farright_jj = farright[jj];
  1237.  
  1238.                 for (c=0; c<channels; c++)
  1239.                   *output_row_ptr++ =
  1240.                     F_ROUND_AND_CLAMP_0255(       coef_left[c] * left_jj
  1241.                                             +   coef_center[c] * center_jj
  1242.                                             +    coef_right[c] * right_jj
  1243.                                             + coef_farright[c] * farright_jj);
  1244.               }
  1245.           }
  1246.  
  1247.         /*
  1248.          * Now, we deal with j = n-2.
  1249.          */
  1250.         for (c=0; c<channels; c++)
  1251.           {
  1252.             a_middle_left[c]      = a_middle_center[c];
  1253.             a_middle_center[c]    = a_middle_right[c];
  1254.             a_middle_right[c]     = a_middle_farright[c];
  1255.  
  1256.             a_bottom_left[c]      = a_bottom_center[c];
  1257.             a_bottom_center[c]    = a_bottom_right[c];
  1258.             a_bottom_right[c]     = a_bottom_farright[c];
  1259.  
  1260.             a_farbottom_left[c]   = a_farbottom_center[c];
  1261.             a_farbottom_center[c] = a_farbottom_right[c];
  1262.             a_farbottom_right[c]  = a_farbottom_farright[c];
  1263.  
  1264.             coef_left[c]   =   middle_ii    * a_middle_left[c]
  1265.                              + bottom_ii    * a_bottom_left[c]
  1266.                              + farbottom_ii * a_farbottom_left[c];
  1267.             coef_center[c] =   middle_ii    * a_middle_center[c]
  1268.                              + bottom_ii    * a_bottom_center[c]
  1269.                              + farbottom_ii * a_farbottom_center[c];
  1270.             coef_right[c]  =   middle_ii    * a_middle_right[c]
  1271.                              + bottom_ii    * a_bottom_right[c]
  1272.                              + farbottom_ii * a_farbottom_right[c];
  1273.           }
  1274.  
  1275.         last_jj = nn_minus_1;
  1276.  
  1277.         for (; jj<last_jj; jj++)
  1278.           {
  1279.             left_jj   =   left[jj];
  1280.             center_jj = center[jj];
  1281.             right_jj  =  right[jj];
  1282.  
  1283.             for (c=0; c<channels; c++)
  1284.               *output_row_ptr++ =
  1285.                 F_ROUND_AND_CLAMP_0255(     coef_left[c] * left_jj
  1286.                                         + coef_center[c] * center_jj
  1287.                                         +  coef_right[c] * right_jj );
  1288.           }
  1289.  
  1290.         /*
  1291.          * Now, we deal with j = n-1 (only one fine pixel left):
  1292.          */
  1293.         for (c=0; c<channels; c++)
  1294.           {
  1295.             a_middle_left[c]      = a_middle_center[c];
  1296.             a_middle_center[c]    = a_middle_right[c];
  1297.  
  1298.             a_bottom_left[c]      = a_bottom_center[c];
  1299.             a_bottom_center[c]    = a_bottom_right[c];
  1300.  
  1301.             a_farbottom_left[c]   = a_farbottom_center[c];
  1302.             a_farbottom_center[c] = a_farbottom_right[c];
  1303.  
  1304.             coef_left[c]   =   middle_ii    * a_middle_left[c]
  1305.                              + bottom_ii    * a_bottom_left[c]
  1306.                              + farbottom_ii * a_farbottom_left[c];
  1307.             coef_center[c] =   middle_ii    * a_middle_center[c]
  1308.                              + bottom_ii    * a_bottom_center[c]
  1309.                              + farbottom_ii * a_farbottom_center[c];
  1310.           }
  1311.  
  1312.         left_jj   =   left[nn_minus_1];
  1313.         center_jj = center[nn_minus_1];
  1314.  
  1315.         for (c=0; c<channels; c++)
  1316.           *output_row_ptr++ =
  1317.             F_ROUND_AND_CLAMP_0255(   coef_left[c]   * left_jj
  1318.                                     + coef_center[c] * center_jj );
  1319.         /*
  1320.          * Write out the row:
  1321.          */
  1322.         gimp_pixel_rgn_set_row (&output_image, output_row, 0, ii, nn);
  1323.       }
  1324.  
  1325.     gimp_progress_update(.005+((gdouble)m*n)/((gdouble)mm*nn)*.17);
  1326.  
  1327.     /*
  1328.      * Now, we deal with the second and later coarse rows, up to the
  1329.      * third to last, inclusive.
  1330.      */
  1331.     for (i=1; i<m_minus_2; i++)
  1332.       {
  1333.         a_temp      = a_top;
  1334.         a_top       = a_middle;
  1335.         a_middle    = a_bottom;
  1336.         a_bottom    = a_farbottom;
  1337.         a_farbottom = a_temp;
  1338.  
  1339.         /*
  1340.          * Row solve on a_farbottom:
  1341.          */
  1342.  
  1343.         a_row_p1 = a_farbottom;
  1344.  
  1345.         for (kk=n_times_channels; kk; kk--, a_row_p1++, a_ptr++)
  1346.           *a_row_p1 = ( (gint) *a_ptr + UNSCALED_SHIFT ) * UNSCALING;
  1347.  
  1348.         /*
  1349.          * Forward substitution:
  1350.          */
  1351.  
  1352.         /*
  1353.          * We'll take care of the very first pixel on the way back.
  1354.          */
  1355.         a_row_p1 = a_farbottom+channels;
  1356.         a_row_p2 = a_farbottom;
  1357.  
  1358.         for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1359.           *a_row_p1 += *a_row_p2 * MINUS_C0_F;
  1360.  
  1361.         for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1362.           *a_row_p1 += *a_row_p2 * MINUS_C1_F;
  1363.  
  1364.         for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1365.           *a_row_p1 += *a_row_p2 * MINUS_C2_F;
  1366.  
  1367.         for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1368.           *a_row_p1 += *a_row_p2 * MINUS_C3_F;
  1369.  
  1370.         for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1371.           *a_row_p1 += *a_row_p2 * MINUS_C4_F;
  1372.  
  1373.         for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1374.           *a_row_p1 += *a_row_p2 * MINUS_C5_F;
  1375.  
  1376.         /*
  1377.          * Because integer overflow is not a risk any more---we do clamp
  1378.          * values later on, and we'd have to even in exact arithmetic---we
  1379.          * use the rounded value of CINFTY instead of its truncated one.
  1380.          */
  1381.         for (kk=n_minus_8_times_channels; kk; kk--, a_row_p1++, a_row_p2++)
  1382.           *a_row_p1 += *a_row_p2 * MINUS_CINFTY_F;
  1383.  
  1384.         for (c=channels; c; c--, a_row_p1++, a_row_p2++)
  1385.           *a_row_p1 = ( *a_row_p2 * MINUS_CINFTY_F + *a_row_p1 ) * CLAST_F;
  1386.  
  1387.         a_row_p1--;
  1388.         a_row_p2 = a_row_p1;
  1389.         for (c=channels; c; c--)
  1390.           a_row_p1--;
  1391.  
  1392.         for (kk=n_minus_7_times_channels; kk; kk--, a_row_p1--, a_row_p2--)
  1393.           *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * CINFTY_F;
  1394.  
  1395.         for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1396.           *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C5_F;
  1397.  
  1398.         for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1399.           *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C4_F;
  1400.  
  1401.         for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1402.           *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C3_F;
  1403.  
  1404.         for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1405.           *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C2_F;
  1406.  
  1407.         for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1408.           *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C1_F;
  1409.  
  1410.         for (c=channels; c; c--, a_row_p1--, a_row_p2--)
  1411.           *a_row_p1 = ( *a_row_p1 - *a_row_p2 ) * C0_F;
  1412.  
  1413.         last_ii = first_past_ii[i];
  1414.  
  1415.         for (; ii<last_ii; ii++)
  1416.           {
  1417.             output_row_ptr = output_row;
  1418.  
  1419.             top_ii       =       top[ii];
  1420.             middle_ii    =    middle[ii];
  1421.             bottom_ii    =    bottom[ii];
  1422.             farbottom_ii = farbottom[ii];
  1423.  
  1424.             a_top_p       = a_top;
  1425.             a_middle_p    = a_middle;
  1426.             a_bottom_p    = a_bottom;
  1427.             a_farbottom_p = a_farbottom;
  1428.  
  1429.             for (c=0; c<channels; c++)
  1430.               {
  1431.                 a_top_center[c]       = *a_top_p++;
  1432.                 a_middle_center[c]    = *a_middle_p++;
  1433.                 a_bottom_center[c]    = *a_bottom_p++;
  1434.                 a_farbottom_center[c] = *a_farbottom_p++;
  1435.               }
  1436.  
  1437.             for (c=0; c<channels; c++)
  1438.               {
  1439.                 a_top_right[c]       = *a_top_p++;
  1440.                 a_middle_right[c]    = *a_middle_p++;
  1441.                 a_bottom_right[c]    = *a_bottom_p++;
  1442.                 a_farbottom_right[c] = *a_farbottom_p++;
  1443.               }
  1444.  
  1445.             for (c=0; c<channels; c++)
  1446.               {
  1447.                 a_top_farright[c]       = *a_top_p++;
  1448.                 a_middle_farright[c]    = *a_middle_p++;
  1449.                 a_bottom_farright[c]    = *a_bottom_p++;
  1450.                 a_farbottom_farright[c] = *a_farbottom_p++;
  1451.               }
  1452.  
  1453.             for (c=0; c<channels; c++)
  1454.               {
  1455.                 coef_center[c]   =         top_ii * a_top_center[c]
  1456.                                    +    middle_ii * a_middle_center[c]
  1457.                                    +    bottom_ii * a_bottom_center[c]
  1458.                                    + farbottom_ii * a_farbottom_center[c];
  1459.                 coef_right[c]    =         top_ii * a_top_right[c]
  1460.                                    +    middle_ii * a_middle_right[c]
  1461.                                    +    bottom_ii * a_bottom_right[c]
  1462.                                    + farbottom_ii * a_farbottom_right[c];
  1463.                 coef_farright[c] =         top_ii * a_top_farright[c]
  1464.                                    +    middle_ii * a_middle_farright[c]
  1465.                                    +    bottom_ii * a_bottom_farright[c]
  1466.                                    + farbottom_ii * a_farbottom_farright[c];
  1467.               }
  1468.  
  1469.             last_jj = first_past_jj[0];
  1470.  
  1471.             for (jj=0; jj<last_jj; jj++)
  1472.               {
  1473.                 center_jj   =   center[jj];
  1474.                 right_jj    =    right[jj];
  1475.                 farright_jj = farright[jj];
  1476.  
  1477.                 for (c=0; c<channels; c++)
  1478.                   *output_row_ptr++ =
  1479.                     F_ROUND_AND_CLAMP_0255(     coef_center[c] * center_jj
  1480.                                             +    coef_right[c] * right_jj
  1481.                                             + coef_farright[c] * farright_jj);
  1482.               }
  1483.  
  1484.             for (j=1; j<n_minus_2; j++)
  1485.               {
  1486.                 for (c=0; c<channels; c++)
  1487.                   {
  1488.                     a_top_left[c]           = a_top_center[c];
  1489.                     a_top_center[c]         = a_top_right[c];
  1490.                     a_top_right[c]          = a_top_farright[c];
  1491.                     a_top_farright[c]       = *a_top_p++;
  1492.  
  1493.                     a_middle_left[c]        = a_middle_center[c];
  1494.                     a_middle_center[c]      = a_middle_right[c];
  1495.                     a_middle_right[c]       = a_middle_farright[c];
  1496.                     a_middle_farright[c]    = *a_middle_p++;
  1497.                    
  1498.                     a_bottom_left[c]        = a_bottom_center[c];
  1499.                     a_bottom_center[c]      = a_bottom_right[c];
  1500.                     a_bottom_right[c]       = a_bottom_farright[c];
  1501.                     a_bottom_farright[c]    = *a_bottom_p++;
  1502.  
  1503.                     a_farbottom_left[c]     = a_farbottom_center[c];
  1504.                     a_farbottom_center[c]   = a_farbottom_right[c];
  1505.                     a_farbottom_right[c]    = a_farbottom_farright[c];
  1506.                     a_farbottom_farright[c] = *a_farbottom_p++;
  1507.  
  1508.                     coef_left[c]     =      top_ii    * a_top_left[c]
  1509.                                        + middle_ii    * a_middle_left[c]
  1510.                                        + bottom_ii    * a_bottom_left[c]
  1511.                                        + farbottom_ii * a_farbottom_left[c];
  1512.                     coef_center[c]   =      top_ii    * a_top_center[c]
  1513.                                        + middle_ii    * a_middle_center[c]
  1514.                                        + bottom_ii    * a_bottom_center[c]
  1515.                                        + farbottom_ii * a_farbottom_center[c];
  1516.                     coef_right[c]    =      top_ii    * a_top_right[c]
  1517.                                        + middle_ii    * a_middle_right[c]
  1518.                                        + bottom_ii    * a_bottom_right[c]
  1519.                                        + farbottom_ii * a_farbottom_right[c];
  1520.                     coef_farright[c] =      top_ii    * a_top_farright[c]
  1521.                                        + middle_ii    * a_middle_farright[c]
  1522.                                        + bottom_ii    * a_bottom_farright[c]
  1523.                                       + farbottom_ii * a_farbottom_farright[c];
  1524.                   }
  1525.  
  1526.                 last_jj = first_past_jj[j];
  1527.  
  1528.                 for (; jj<last_jj; jj++)
  1529.                   {
  1530.                     left_jj     =     left[jj];
  1531.                     center_jj   =   center[jj];
  1532.                     right_jj    =    right[jj];
  1533.                     farright_jj = farright[jj];
  1534.  
  1535.                     for (c=0; c<channels; c++)
  1536.                       *output_row_ptr++ =
  1537.                         F_ROUND_AND_CLAMP_0255(       coef_left[c] * left_jj
  1538.                                                 +   coef_center[c] * center_jj
  1539.                                                 +    coef_right[c] * right_jj
  1540.                                             + coef_farright[c] * farright_jj );
  1541.                   }
  1542.               }
  1543.  
  1544.             /*
  1545.              * Now, we deal with j = n-2.
  1546.              */
  1547.             for (c=0; c<channels; c++)
  1548.               {
  1549.                 a_top_left[c]         = a_top_center[c];
  1550.                 a_top_center[c]       = a_top_right[c];
  1551.                 a_top_right[c]        = a_top_farright[c];
  1552.                
  1553.                 a_middle_left[c]      = a_middle_center[c];
  1554.                 a_middle_center[c]    = a_middle_right[c];
  1555.                 a_middle_right[c]     = a_middle_farright[c];
  1556.  
  1557.                 a_bottom_left[c]      = a_bottom_center[c];
  1558.                 a_bottom_center[c]    = a_bottom_right[c];
  1559.                 a_bottom_right[c]     = a_bottom_farright[c];
  1560.  
  1561.                 a_farbottom_left[c]   = a_farbottom_center[c];
  1562.                 a_farbottom_center[c] = a_farbottom_right[c];
  1563.                 a_farbottom_right[c]  = a_farbottom_farright[c];
  1564.  
  1565.                 coef_left[c]   =         top_ii * a_top_left[c]
  1566.                                  +    middle_ii * a_middle_left[c]
  1567.                                  +    bottom_ii * a_bottom_left[c]
  1568.                                  + farbottom_ii * a_farbottom_left[c];
  1569.                 coef_center[c] =         top_ii * a_top_center[c]
  1570.                                  +    middle_ii * a_middle_center[c]
  1571.                                  +    bottom_ii * a_bottom_center[c]
  1572.                                  + farbottom_ii * a_farbottom_center[c];
  1573.                 coef_right[c]  =         top_ii * a_top_right[c]
  1574.                                  +    middle_ii * a_middle_right[c]
  1575.                                  +    bottom_ii * a_bottom_right[c]
  1576.                                  + farbottom_ii * a_farbottom_right[c];
  1577.               }
  1578.  
  1579.             last_jj = nn_minus_1;
  1580.        
  1581.             for (; jj<last_jj; jj++)
  1582.               {
  1583.                 left_jj   =   left[jj];
  1584.                 center_jj = center[jj];
  1585.                 right_jj  =  right[jj];
  1586.  
  1587.                 for (c=0; c<channels; c++)
  1588.                   *output_row_ptr++ =
  1589.                     F_ROUND_AND_CLAMP_0255(     coef_left[c] * left_jj
  1590.                                             + coef_center[c] * center_jj
  1591.                                             +  coef_right[c] * right_jj );
  1592.               }
  1593.  
  1594.             /*
  1595.              * Now, we deal with j = n-1 (the very last fine pixel).
  1596.              */
  1597.             for (c=0; c<channels; c++)
  1598.               {
  1599.                 a_top_left[c]         = a_top_center[c];
  1600.                 a_top_center[c]       = a_top_right[c];
  1601.  
  1602.                 a_middle_left[c]      = a_middle_center[c];
  1603.                 a_middle_center[c]    = a_middle_right[c];
  1604.            
  1605.                 a_bottom_left[c]      = a_bottom_center[c];
  1606.                 a_bottom_center[c]    = a_bottom_right[c];
  1607.  
  1608.                 a_farbottom_left[c]   = a_farbottom_center[c];
  1609.                 a_farbottom_center[c] = a_farbottom_right[c];
  1610.  
  1611.                 coef_left[c]   =         top_ii * a_top_left[c]
  1612.                                  +    middle_ii * a_middle_left[c]
  1613.                                  +    bottom_ii * a_bottom_left[c]
  1614.                                  + farbottom_ii * a_farbottom_left[c];
  1615.                 coef_center[c] =         top_ii * a_top_center[c]
  1616.                                  +    middle_ii * a_middle_center[c]
  1617.                                  +    bottom_ii * a_bottom_center[c]
  1618.                                  + farbottom_ii * a_farbottom_center[c];
  1619.               }
  1620.  
  1621.             left_jj   =   left[nn_minus_1];
  1622.             center_jj = center[nn_minus_1];
  1623.            
  1624.             for (c=0; c<channels; c++)
  1625.               *output_row_ptr++ =
  1626.                 F_ROUND_AND_CLAMP_0255(     coef_left[c] * left_jj
  1627.                                         + coef_center[c] * center_jj );
  1628.  
  1629.             gimp_pixel_rgn_set_row (&output_image, output_row, 0, ii, nn);
  1630.           }
  1631.       }
  1632.  
  1633.     gimp_progress_update(.99);
  1634.  
  1635.     /*
  1636.      * Now, we deal with the second to last coarse row:
  1637.      */
  1638.     a_temp      = a_top;
  1639.     a_top       = a_middle;
  1640.     a_middle    = a_bottom;
  1641.     a_bottom    = a_farbottom;
  1642.     a_farbottom = a_temp;
  1643.  
  1644.     last_ii = mm - 1;
  1645.  
  1646.     for (; ii<last_ii; ii++)
  1647.       {
  1648.         output_row_ptr = output_row;
  1649.  
  1650.         top_ii    =    top[ii];
  1651.         middle_ii = middle[ii];
  1652.         bottom_ii = bottom[ii];
  1653.  
  1654.         a_top_p    = a_top;
  1655.         a_middle_p = a_middle;
  1656.         a_bottom_p = a_bottom;
  1657.  
  1658.         for (c=0; c<channels; c++)
  1659.           {
  1660.             a_top_center[c]      = *a_top_p++;
  1661.             a_middle_center[c]   = *a_middle_p++;
  1662.             a_bottom_center[c]   = *a_bottom_p++;
  1663.           }
  1664.  
  1665.         for (c=0; c<channels; c++)
  1666.           {
  1667.             a_top_right[c]       = *a_top_p++;
  1668.             a_middle_right[c]    = *a_middle_p++;
  1669.             a_bottom_right[c]    = *a_bottom_p++;
  1670.           }
  1671.  
  1672.         for (c=0; c<channels; c++)
  1673.           {
  1674.             a_top_farright[c]    = *a_top_p++;
  1675.             a_middle_farright[c] = *a_middle_p++;
  1676.             a_bottom_farright[c] = *a_bottom_p++;
  1677.           }
  1678.  
  1679.         for (c=0; c<channels; c++)
  1680.           {
  1681.             coef_center[c]   =      top_ii * a_top_center[c]
  1682.                                + middle_ii * a_middle_center[c]
  1683.                                + bottom_ii * a_bottom_center[c];
  1684.             coef_right[c]    =      top_ii * a_top_right[c]
  1685.                                + middle_ii * a_middle_right[c]
  1686.                                + bottom_ii * a_bottom_right[c];
  1687.             coef_farright[c] =      top_ii * a_top_farright[c]
  1688.                                + middle_ii * a_middle_farright[c]
  1689.                                + bottom_ii * a_bottom_farright[c];
  1690.           }
  1691.        
  1692.         last_jj = first_past_jj[0];
  1693.  
  1694.         for (jj=0; jj<last_jj; jj++)
  1695.           {
  1696.             center_jj   =   center[jj];
  1697.             right_jj    =    right[jj];
  1698.             farright_jj = farright[jj];
  1699.  
  1700.             for (c=0; c<channels; c++)
  1701.               *output_row_ptr++ =
  1702.                 F_ROUND_AND_CLAMP_0255(   coef_center[c]   * center_jj
  1703.                                          + coef_right[c]    * right_jj
  1704.                                       + coef_farright[c] * farright_jj );
  1705.           }
  1706.  
  1707.         for (j=1; j<n_minus_2; j++)
  1708.           {
  1709.             for (c=0; c<channels; c++)
  1710.               {
  1711.                 a_top_left[c]        = a_top_center[c];
  1712.                 a_top_center[c]      = a_top_right[c];
  1713.                 a_top_right[c]       = a_top_farright[c];
  1714.                 a_top_farright[c]    = *a_top_p++;
  1715.  
  1716.                 a_middle_left[c]     = a_middle_center[c];
  1717.                 a_middle_center[c]   = a_middle_right[c];
  1718.                 a_middle_right[c]    = a_middle_farright[c];
  1719.                 a_middle_farright[c] = *a_middle_p++;
  1720.  
  1721.                 a_bottom_left[c]     = a_bottom_center[c];
  1722.                 a_bottom_center[c]   = a_bottom_right[c];
  1723.                 a_bottom_right[c]    = a_bottom_farright[c];
  1724.                 a_bottom_farright[c] = *a_bottom_p++;
  1725.  
  1726.                 coef_left[c]     =      top_ii * a_top_left[c]
  1727.                                    + middle_ii * a_middle_left[c]
  1728.                                    + bottom_ii * a_bottom_left[c];
  1729.                 coef_center[c]   =      top_ii * a_top_center[c]
  1730.                                    + middle_ii * a_middle_center[c]
  1731.                                    + bottom_ii * a_bottom_center[c];
  1732.                 coef_right[c]    =      top_ii * a_top_right[c]
  1733.                                    + middle_ii * a_middle_right[c]
  1734.                                    + bottom_ii * a_bottom_right[c];
  1735.                 coef_farright[c] =      top_ii * a_top_farright[c]
  1736.                                    + middle_ii * a_middle_farright[c]
  1737.                                    + bottom_ii * a_bottom_farright[c];
  1738.               }
  1739.  
  1740.             last_jj = first_past_jj[j];
  1741.            
  1742.             for (; jj<last_jj; jj++)
  1743.               {
  1744.                 left_jj     =     left[jj];
  1745.                 center_jj   =   center[jj];
  1746.                 right_jj    =    right[jj];
  1747.                 farright_jj = farright[jj];
  1748.                
  1749.                 for (c=0; c<channels; c++)
  1750.                   *output_row_ptr++ =
  1751.                     F_ROUND_AND_CLAMP_0255(       coef_left[c] * left_jj
  1752.                                             +   coef_center[c] * center_jj
  1753.                                             +    coef_right[c] * right_jj
  1754.                                             + coef_farright[c] * farright_jj );
  1755.               }
  1756.           }
  1757.  
  1758.         /*
  1759.          * Now, we deal with j = n-2:
  1760.          */
  1761.         for (c=0; c<channels; c++)
  1762.           {
  1763.             a_top_left[c]      = a_top_center[c];
  1764.             a_top_center[c]    = a_top_right[c];
  1765.             a_top_right[c]     = a_top_farright[c];
  1766.            
  1767.             a_middle_left[c]   = a_middle_center[c];
  1768.             a_middle_center[c] = a_middle_right[c];
  1769.             a_middle_right[c]  = a_middle_farright[c];
  1770.  
  1771.             a_bottom_left[c]   = a_bottom_center[c];
  1772.             a_bottom_center[c] = a_bottom_right[c];
  1773.             a_bottom_right[c]  = a_bottom_farright[c];
  1774.  
  1775.             coef_left[c]   =      top_ii * a_top_left[c]
  1776.                              + middle_ii * a_middle_left[c]
  1777.                              + bottom_ii * a_bottom_left[c];
  1778.             coef_center[c] =      top_ii * a_top_center[c]
  1779.                              + middle_ii * a_middle_center[c]
  1780.                              + bottom_ii * a_bottom_center[c];
  1781.             coef_right[c]  =      top_ii * a_top_right[c]
  1782.                              + middle_ii * a_middle_right[c]
  1783.                              + bottom_ii * a_bottom_right[c];
  1784.           }
  1785.  
  1786.         last_jj = nn_minus_1;
  1787.  
  1788.         for (; jj<last_jj; jj++)
  1789.           {
  1790.             left_jj   =   left[jj];
  1791.             center_jj = center[jj];
  1792.             right_jj  =  right[jj];
  1793.  
  1794.             for (c=0; c<channels; c++)
  1795.               *output_row_ptr++ =
  1796.                 F_ROUND_AND_CLAMP_0255(     coef_left[c] * left_jj
  1797.                                         + coef_center[c] * center_jj
  1798.                                         +  coef_right[c] * right_jj );
  1799.           }
  1800.  
  1801.         /*
  1802.          * Now, we deal with j = n-1 (last fine pixel):
  1803.          */
  1804.         for (c=0; c<channels; c++)
  1805.           {
  1806.             a_top_left[c]      = a_top_center[c];
  1807.             a_top_center[c]    = a_top_right[c];
  1808.            
  1809.             a_middle_left[c]   = a_middle_center[c];
  1810.             a_middle_center[c] = a_middle_right[c];
  1811.            
  1812.             a_bottom_left[c]   = a_bottom_center[c];
  1813.             a_bottom_center[c] = a_bottom_right[c];
  1814.            
  1815.             coef_left[c]   =      top_ii * a_top_left[c]
  1816.                              + middle_ii * a_middle_left[c]
  1817.                              + bottom_ii * a_bottom_left[c];
  1818.             coef_center[c] =      top_ii * a_top_center[c]
  1819.                              + middle_ii * a_middle_center[c]
  1820.                              + bottom_ii * a_bottom_center[c];
  1821.           }
  1822.  
  1823.         left_jj   =   left[nn_minus_1];
  1824.         center_jj = center[nn_minus_1];
  1825.  
  1826.         for (c=0; c<channels; c++)
  1827.           *output_row_ptr++ =
  1828.             F_ROUND_AND_CLAMP_0255(     coef_left[c] * left_jj
  1829.                                     + coef_center[c] * center_jj );
  1830.  
  1831.         gimp_pixel_rgn_set_row (&output_image, output_row, 0, ii, nn);
  1832.       }
  1833.  
  1834.     /*
  1835.      * Compute the last coarse row:
  1836.      */
  1837.     a_temp   = a_top;
  1838.     a_top    = a_middle;
  1839.     a_middle = a_bottom;
  1840.     a_bottom = a_temp;
  1841.    
  1842.     output_row_ptr = output_row;
  1843.  
  1844.     top_ii    = top[ii];
  1845.     middle_ii = middle[ii];
  1846.  
  1847.     a_top_p    = a_top;
  1848.     a_middle_p = a_middle;
  1849.  
  1850.     for (c=0; c<channels; c++)
  1851.       {
  1852.         a_top_center[c]      = *a_top_p++;
  1853.         a_middle_center[c]   = *a_middle_p++;
  1854.       }
  1855.  
  1856.     for (c=0; c<channels; c++)
  1857.       {
  1858.         a_top_right[c]       = *a_top_p++;
  1859.         a_middle_right[c]    = *a_middle_p++;
  1860.       }
  1861.  
  1862.     for (c=0; c<channels; c++)
  1863.       {
  1864.         a_top_farright[c]    = *a_top_p++;
  1865.         a_middle_farright[c] = *a_middle_p++;
  1866.       }
  1867.  
  1868.     for (c=0; c<channels; c++)
  1869.       {
  1870.         coef_center[c]   =      top_ii * a_top_center[c]
  1871.                            + middle_ii * a_middle_center[c];
  1872.         coef_right[c]    =      top_ii * a_top_right[c]
  1873.                            + middle_ii * a_middle_right[c];
  1874.         coef_farright[c] =      top_ii * a_top_farright[c]
  1875.                            + middle_ii * a_middle_farright[c];
  1876.       }
  1877.  
  1878.     last_jj = first_past_jj[0];
  1879.  
  1880.     for (jj=0; jj<last_jj; jj++)
  1881.       {
  1882.         center_jj   =   center[jj];
  1883.         right_jj    =    right[jj];
  1884.         farright_jj = farright[jj];
  1885.  
  1886.         for (c=0; c<channels; c++)
  1887.           *output_row_ptr++ =
  1888.             F_ROUND_AND_CLAMP_0255(   coef_center[c]   * center_jj
  1889.                                     + coef_right[c]    * right_jj
  1890.                                     + coef_farright[c] * farright_jj );
  1891.       }
  1892.  
  1893.     for (j=1; j<n_minus_2; j++)
  1894.       {
  1895.         for (c=0; c<channels; c++)
  1896.           {
  1897.             a_top_left[c]        = a_top_center[c];
  1898.             a_top_center[c]      = a_top_right[c];
  1899.             a_top_right[c]       = a_top_farright[c];
  1900.             a_top_farright[c]    = *a_top_p++;
  1901.  
  1902.             a_middle_left[c]     = a_middle_center[c];
  1903.             a_middle_center[c]   = a_middle_right[c];
  1904.             a_middle_right[c]    = a_middle_farright[c];
  1905.             a_middle_farright[c] = *a_middle_p++;
  1906.            
  1907.             coef_left[c]     =      top_ii * a_top_left[c]
  1908.                                + middle_ii * a_middle_left[c];
  1909.             coef_center[c]   =      top_ii * a_top_center[c]
  1910.                                + middle_ii * a_middle_center[c];
  1911.             coef_right[c]    =      top_ii * a_top_right[c]
  1912.                                + middle_ii * a_middle_right[c];
  1913.             coef_farright[c] =      top_ii * a_top_farright[c]
  1914.                                + middle_ii * a_middle_farright[c];
  1915.           }
  1916.  
  1917.         last_jj = first_past_jj[j];
  1918.  
  1919.         for (; jj<last_jj; jj++)
  1920.           {
  1921.             left_jj     =     left[jj];
  1922.             center_jj   =   center[jj];
  1923.             right_jj    =    right[jj];
  1924.             farright_jj = farright[jj];
  1925.  
  1926.             for (c=0; c<channels; c++)
  1927.               *output_row_ptr++ =
  1928.                 F_ROUND_AND_CLAMP_0255(       coef_left[c] * left_jj
  1929.                                         +   coef_center[c] * center_jj
  1930.                                         +    coef_right[c] * right_jj
  1931.                                         + coef_farright[c] * farright_jj );
  1932.           }
  1933.       }
  1934.  
  1935.     /*
  1936.      * Now, we deal with j = n-2.
  1937.      */
  1938.     for (c=0; c<channels; c++)
  1939.       {
  1940.         a_top_left[c]      = a_top_center[c];
  1941.         a_top_center[c]    = a_top_right[c];
  1942.         a_top_right[c]     = a_top_farright[c];
  1943.        
  1944.         a_middle_left[c]   = a_middle_center[c];
  1945.         a_middle_center[c] = a_middle_right[c];
  1946.         a_middle_right[c]  = a_middle_farright[c];
  1947.        
  1948.         coef_left[c]   =      top_ii * a_top_left[c]
  1949.                          + middle_ii * a_middle_left[c];
  1950.         coef_center[c] =      top_ii * a_top_center[c]
  1951.                          + middle_ii * a_middle_center[c];
  1952.         coef_right[c]  =      top_ii * a_top_right[c]
  1953.                          + middle_ii * a_middle_right[c];
  1954.       }
  1955.  
  1956.     last_jj = nn_minus_1;
  1957.  
  1958.     for (; jj<last_jj; jj++)
  1959.       {
  1960.         left_jj   =   left[jj];
  1961.         center_jj = center[jj];
  1962.         right_jj  =  right[jj];
  1963.        
  1964.         for (c=0; c<channels; c++)
  1965.           *output_row_ptr++ =
  1966.             F_ROUND_AND_CLAMP_0255(     coef_left[c] * left_jj
  1967.                                     + coef_center[c] * center_jj
  1968.                                     +  coef_right[c] * right_jj );
  1969.       }
  1970.  
  1971.     /*
  1972.      * Now, we deal with j = n-1:
  1973.      */
  1974.     for (c=0; c<channels; c++)
  1975.       {
  1976.         a_top_left[c]      = a_top_center[c];
  1977.         a_top_center[c]    = a_top_right[c];
  1978.        
  1979.         a_middle_left[c]   = a_middle_center[c];
  1980.         a_middle_center[c] = a_middle_right[c];
  1981.          
  1982.         coef_left[c]   =      top_ii * a_top_left[c]
  1983.           + middle_ii * a_middle_left[c];
  1984.         coef_center[c] =      top_ii * a_top_center[c]
  1985.           + middle_ii * a_middle_center[c];
  1986.       }
  1987.  
  1988.     left_jj   =   left[nn_minus_1];
  1989.     center_jj = center[nn_minus_1];
  1990.  
  1991.     for (c=0; c<channels; c++)
  1992.       *output_row_ptr++ =
  1993.         F_ROUND_AND_CLAMP_0255(     coef_left[c] * left_jj
  1994.                                 + coef_center[c] * center_jj );
  1995.  
  1996.     gimp_pixel_rgn_set_row (&output_image, output_row, 0, ii, nn);
  1997.  
  1998.     g_free(a_farbottom);
  1999.     g_free(a_bottom);
  2000.     g_free(a_middle);
  2001.     g_free(a_top);
  2002.     g_free(output_row);
  2003.   }
  2004.  
  2005.   g_free(a);
  2006.  
  2007.   gimp_drawable_flush (new_drawable);
  2008.   gimp_drawable_merge_shadow (new_drawable->drawable_id, FALSE);
  2009.   gimp_drawable_update (new_drawable->drawable_id, 0, 0, nn, mm);
  2010.   gimp_drawable_detach (new_drawable);
  2011.  
  2012.   return new_image_ID;
  2013. }
  2014.  
  2015. static void
  2016. first_past_index (gint o,
  2017.                   gint oo,
  2018.                   gint first_past_kk[])
  2019. {
  2020.   /*
  2021.    * The comments use the variables which are relevant to the use of
  2022.    * this function in the horizontal direction. That is, n is o, nn is
  2023.    * oo, j is k, and jj is kk.
  2024.    *
  2025.    * This subprogram puts in first_past_jj[j], for j in 0..n-2, the
  2026.    * index of the first fine cell with a "center" which is not to the
  2027.    * left of the center of the j+1th coarse cell (the last cell has
  2028.    * index n-1, and so has no right neighbour). An alternate
  2029.    * description: first_past_jj[j] is the index of the first fine cell
  2030.    * with center >= j+3/2. This index matters for the following
  2031.    * reason: It is the index at which the rightmost relevant B-Spline
  2032.    * coefficient shifts to the right (center = j+3/2 could go either
  2033.    * way).
  2034.    *
  2035.    * This code uses a general convention that the input pixels have
  2036.    * unit sides, and that (describing things in 1D) the global
  2037.    * coordinate of the left endpoint of an input cell is exactly its
  2038.    * index. That is, the first coarse (input) cell goes from 0 to 1,
  2039.    * the last coarse cell from n-1 to n.
  2040.    *
  2041.    * Because of this convention, there is an alternate description of
  2042.    * first_past_jj[j]: Because coarse cell centers are 1/2 to the
  2043.    * right of left boundaries, first_past_jj is also the index of the
  2044.    * first fine cell for which the left end of the box filtering mask
  2045.    * (which has diameter 1/2) is >= j+1.  The very first "box" starts
  2046.    * at zero, and it shifts by h at every index shift. This gives jj*h
  2047.    * >= j+1, that is, it is the first jj which violates jj*(n-1) <
  2048.    * (j+1)*(nn-1).
  2049.    *
  2050.    * Because the jth cell is [j,j+1], and choosing the very first fine
  2051.    * cell center to be at 1/2 (same as the very first coarse cell
  2052.    * center), so that the very last fine, and coarse, cell center is
  2053.    * n-1/2, we have that as the fine index runs from 0 to nn-1, a
  2054.    * distance of n-1 is covered, so that
  2055.    *
  2056.    * h = (n-1)/(nn-1)
  2057.    *
  2058.    * (the usual step for interpolation methods). Consequently,
  2059.    * first_past_jj[j] is the first index jj for which
  2060.    *
  2061.    * 1/2 + jj * h >= j + 3/2;
  2062.    *
  2063.    * in integer arithmetic:
  2064.    *
  2065.    * jj * (n-1) >= (j+1) * (nn-1).
  2066.    *
  2067.    * Alternately, it is the first index for which the condition
  2068.    *
  2069.    * jj * (n-1) < (j+1) * (nn-1)
  2070.    *
  2071.    * does not hold.
  2072.    */
  2073.  
  2074.   gint o_minus_1 = o-1;
  2075.  
  2076.   gint k = 0;
  2077.  
  2078.   if (oo>o)
  2079.   {
  2080.     gint kk = 1;
  2081.     gint kk_times_o_minus_1 = o_minus_1;
  2082.     gint oo_minus_1 = oo-1;
  2083.     gint k_plus_1_times_oo_minus_1 = oo_minus_1;
  2084.  
  2085.     do {
  2086.       while (kk_times_o_minus_1 < k_plus_1_times_oo_minus_1)
  2087.         {
  2088.           kk++;
  2089.           kk_times_o_minus_1 += o_minus_1;
  2090.         }
  2091.       first_past_kk[k]=kk;
  2092.       kk++;
  2093.       kk_times_o_minus_1 += o_minus_1;
  2094.       k++;
  2095.       k_plus_1_times_oo_minus_1 += oo_minus_1;
  2096.     } while (k<o_minus_1);
  2097.   }
  2098.   else
  2099.   {
  2100.     do
  2101.       {
  2102.         first_past_kk[k] = k+1;
  2103.       }
  2104.     while (++k<o_minus_1);
  2105.   }
  2106. }
  2107.  
  2108. /*
  2109.  * smooth_coarse_to_fine_coefficients ()
  2110.  *
  2111.  * This function is called twice: once in the "vertical" direction,
  2112.  * and once in the "horizontal" direction.
  2113.  *
  2114.  * The comments use the variables which are relevant to the use of
  2115.  * this function in the horizontal direction. That is, n is o, nn is
  2116.  * oo, j is k, and jj is kk.
  2117.  *
  2118.  * "INPUTS:" number of coarse cells n
  2119.  *           number of fine cells nn
  2120.  *           one int array with at least n entries: first_past_jj[]
  2121.  *
  2122.  * "OUTPUTS:" four double arrays with at least nn entries:
  2123.  *            first[], second[], third[], fourth[]
  2124.  */
  2125. static void
  2126. smooth_coarse_to_fine_coefficients (gint   o,
  2127.                                     gint   oo,
  2128.                                     gint   first_past_kk[],
  2129.                                     gfloat first[],
  2130.                                     gfloat second[],
  2131.                                     gfloat third[],
  2132.                                     gfloat fourth[])
  2133. {
  2134.   /*
  2135.    * Possible improvements: use long or unsigned integers to prevent
  2136.    * integer overflow: integers as large as n times nn are computed.
  2137.    * This does not appear to have been a problem, but it could be.
  2138.    */
  2139.  
  2140.   /*
  2141.    * Possible improvement: Use left/right symmetry to eliminate the
  2142.    * computation of half the coefficients. This is a Theta(m+n)
  2143.    * improvement (probably not worth the trouble). Because of the way
  2144.    * indices "stay left" for fine cells which overlap two coarse
  2145.    * cells, this is not as simple as it looks, unless there are no
  2146.    * overlapping fine cells, that is, unless the magnification factor
  2147.    * is an integer.
  2148.    */
  2149.  
  2150.   /*
  2151.    * Formulas for the antiderivatives of the cardinal basis functions,
  2152.    * with constants of integration set so that they are equal to zero
  2153.    * at the left boundary of the cell of interest: They are structured
  2154.    * a la Horner to minimize flops. Because x is Theta(1), some mild
  2155.    * cancellation will occur.
  2156.    */
  2157.  
  2158.   /*
  2159.    * The following describe the centered B-Spline (applicable in the
  2160.    * interior):
  2161.    */
  2162. #define LEFT_BSPLINE(x)            ( (x) *         (x) *       (x))
  2163. #define CENTER_BSPLINE(x)          ( (x) * (  3. - (x) * (-3.+(x)+(x))))
  2164. #define RIGHT_BSPLINE(x)           ( (x) * (  3. + (x) * (-3.+(x))))
  2165.   /*
  2166.    * The coordinate x is with respect to the left boundary point of
  2167.    * the applicable cell. That is: x = 1 for LEFT_BSPLINE corresponds
  2168.    * to x = 0 for CENTER_BSPLINE, and x = 1 for CENTER_BSPLINE
  2169.    * corresponds to x = 0 for RIGHT_BSPLINE.
  2170.    */
  2171.  
  2172.   /*
  2173.    * For natural boundary conditions, LEFT_BDRY_LEFT_SPLINE and
  2174.    * RIGHT_BDRY_RIGHT_SPLINE are identical to those which apply in the
  2175.    * interior. Consequently, they are redundant and could be removed
  2176.    * from the code.
  2177.    */
  2178. #define LEFT_BDRY_SPLINE(x)        ( (x) * (  6. -         (x)*(x)))
  2179. #define LEFT_BDRY_LEFT_SPLINE(x)   ( (x) *                 (x)*(x))
  2180. #define RIGHT_BDRY_SPLINE(x)       ( (x) * (  3. + (x) * ( 3.- (x))))
  2181. #define RIGHT_BDRY_RIGHT_SPLINE(x) ( (x) * (  3. + (x) * (-3.+ (x))))
  2182.  
  2183.  /*
  2184.   * jj * (n-1) / (nn-1) (in rational arithmetic)
  2185.   *
  2186.   * is the absolute coordinate of the left endpoint of the index jj
  2187.   * fine cell (recall that jj starts at zero).
  2188.   *
  2189.   * The distance between fine cell centers (and fine cell integration
  2190.   * left boundaries, as well as right boundaries) is (n-1)/(nn-1).
  2191.   *
  2192.   * Now, when the fine cell center is at a coarse cell center or just
  2193.   * past it, the coordinate, with respect to the left boundary of the
  2194.   * coarse cell, of the left boundary point of the integration mask is
  2195.   *
  2196.   * jj * (n-1) / (nn-1) - floor of this value.
  2197.   *
  2198.   * Now, first_past[j-1] is the first jj for which the left boundary
  2199.   * point of the integration mask is not to the left of j.
  2200.   * Consequently, this is the desired floor. This implies that we set
  2201.   *
  2202.   * x = jj * (n-1) / (nn-1) - j
  2203.   *
  2204.   * when we start iterating at the jj value equal to the previous
  2205.   * first_past, which happens to be the value of j at loop entry. In
  2206.   * order to minimize round off, we actually rewrite as
  2207.   *
  2208.   * x = ( jj * (n-1) - j * (nn-1) ) / (nn-1).
  2209.   */
  2210.  
  2211.   /*
  2212.    * For convenience, three coefficients---the last one equal to
  2213.    * zero---are given for the very first pixel in each direction, and
  2214.    * only two for the final pixel in each direction. For the second
  2215.    * and second to last coarse cells, three are given. Otherwise, four
  2216.    * coefficients are given.
  2217.    */
  2218.  
  2219.   gint k;
  2220.   gint kk;
  2221.  
  2222.   gint o_minus_1  = o-1;
  2223.   gint o_minus_2  = o_minus_1 - 1;
  2224.   gint oo_minus_1 = oo-1;
  2225.  
  2226.   gdouble one_over_oo_minus_1 = 1./(oo_minus_1);
  2227.   gdouble h = (o-1) * one_over_oo_minus_1;
  2228.  
  2229.   gdouble x = h;
  2230.  
  2231.   gint first_past = first_past_kk[0];
  2232.  
  2233.   /* first[0] = 0.f; */
  2234.   second[0] = 5.f;
  2235.   third[0]  = 1.f;
  2236.   fourth[0] = 0.f;
  2237.  
  2238.   for (kk=1; kk<first_past; kk++, x+=h)
  2239.     {
  2240.       /*    first[kk]  = 0.f; */
  2241.       /* second[kk] = ( 5. - LEFT_BDRY_SPLINE(x) ) + RIGHT_BSPLINE(x); */
  2242.       second[kk] = 5. + x * ( -3. + x * ( -3. + x + x ) );
  2243.       /* third[kk]  = ( 1. - LEFT_BSPLINE(x) ) + CENTER_BSPLINE(x); */
  2244.       third[kk]  = 1. + 3. * x * ( 1. + x * ( 1. - x ) );
  2245.       fourth[kk] = LEFT_BSPLINE(x);
  2246.     }
  2247.  
  2248.   for (k=1; k<o_minus_2; k++)
  2249.     {
  2250.       x = ( kk*o_minus_1 - k*oo_minus_1 ) * one_over_oo_minus_1;
  2251.       first_past = first_past_kk[k];
  2252.      
  2253.       for (; kk<first_past; kk++, x+=h)
  2254.       {
  2255.         /* first[kk]  = 1. - RIGHT_BSPLINE(x); */
  2256.         first[kk]  = 1. + x * ( -3. + x * ( 3. - x ) );
  2257.         /* second[kk] = ( 4. - CENTER_BSPLINE(x) ) + RIGHT_BSPLINE(x); */
  2258.         second[kk] = 4. + 3. * ( x * x ) * ( -2. + x );
  2259.         /* third[kk]  = ( 1. - LEFT_BSPLINE(x) ) + CENTER_BSPLINE(x); */
  2260.         third[kk]  = 1. + 3. * ( x * ( 1. + x * ( 1. - x ) ) );
  2261.         fourth[kk] = LEFT_BSPLINE(x);
  2262.       }
  2263.     }
  2264.  
  2265.   x = ( kk*o_minus_1 - (o_minus_2)*oo_minus_1 ) * one_over_oo_minus_1;
  2266.   first_past = oo - 1;
  2267.  
  2268.   for (; kk<first_past; kk++, x+=h)
  2269.     {
  2270.       /* first[kk]  = 1. - RIGHT_BSPLINE(x); */
  2271.       first[kk]  = 1. + x * ( -3. + x * ( 3. - x ) );
  2272.       /* second[kk] = ( 4. - CENTER_BSPLINE(x) ) + RIGHT_BSPLINE(x); */
  2273.       second[kk] = 4. + 3. * ( ( x * x ) * ( -2. + x ) );
  2274.       /* third[kk]  = ( 1. - LEFT_BSPLINE(x) ) + RIGHT_BDRY_SPLINE(x); */
  2275.       third[kk]  = 1. + x * ( 3. + x * ( 3. - ( x + x ) ) );
  2276.       /*      fourth[kk] = 0.f; */
  2277.     }
  2278.  
  2279.   first[first_past]  = 1.f;
  2280.   second[first_past] = 5.f;
  2281.   /*  third[first_past]  = 0.f; */
  2282.   /*  fourth[first_past] = 0.f; */
  2283. }
  2284.  
  2285. static gboolean
  2286. scale_dialog (gint32 image_ID,
  2287.               gint   m,
  2288.               gint   n)
  2289. {
  2290.   GtkWidget *dialog;
  2291.   GtkWidget *main_vbox;
  2292.   GtkWidget *main_hbox;
  2293.   GtkWidget *frame;
  2294.   GtkWidget *frame_label;
  2295.   GtkWidget *alignment;
  2296.   GtkWidget *size_entry;
  2297.   GtkWidget *aspect_chain;
  2298.   GimpUnit   unit;
  2299.   gdouble    xres;
  2300.   gdouble    yres;
  2301.   gint       run;
  2302.   gboolean   not_done;
  2303.  
  2304.   gimp_ui_init ("IBFNBQH", FALSE);
  2305.  
  2306.   dialog = gimp_dialog_new ("IBFNBQH", "IBFNBQH",
  2307.                             NULL, 0,
  2308.                             gimp_standard_help_func, "IBFNBQH",
  2309.                             GIMP_STOCK_RESET, RESPONSE_RESET,
  2310.                             GTK_STOCK_CANCEL, GTK_RESPONSE_CANCEL,
  2311.                             GTK_STOCK_OK,     GTK_RESPONSE_OK,
  2312.                             NULL);
  2313.  
  2314.   main_vbox = gtk_vbox_new (FALSE, 12);
  2315.   gtk_container_add (GTK_CONTAINER (GTK_DIALOG (dialog)->vbox), main_vbox);
  2316.   gtk_widget_show (main_vbox);
  2317.  
  2318.   frame = gtk_frame_new (NULL);
  2319.   gtk_widget_show (frame);
  2320.   gtk_box_pack_start (GTK_BOX (main_vbox), frame, TRUE, TRUE, 0);
  2321.   gtk_container_set_border_width (GTK_CONTAINER (frame), 6);
  2322.  
  2323.   alignment = gtk_alignment_new (0.5f, 0.5f, 0.0f, 0.0f);
  2324.   gtk_widget_show (alignment);
  2325.   gtk_container_add (GTK_CONTAINER (frame), alignment);
  2326.   gtk_alignment_set_padding (GTK_ALIGNMENT (alignment), 6, 6, 6, 6);
  2327.  
  2328.   main_hbox = gtk_hbox_new (FALSE, 12);
  2329.   gtk_widget_show (main_hbox);
  2330.   gtk_container_add (GTK_CONTAINER (alignment), main_hbox);
  2331.  
  2332.   /*  Get the image resolution and unit  */
  2333.   gimp_image_get_resolution (image_ID, &xres, &yres);
  2334.   unit = gimp_image_get_unit (image_ID);
  2335.  
  2336.   size_entry = gimp_size_entry_new(2, unit, "%a", TRUE, TRUE, FALSE,
  2337.                                    8,
  2338.                                    GIMP_SIZE_ENTRY_UPDATE_SIZE);
  2339.  
  2340.   /*  set the unit back to pixels, since most times we will want pixels */
  2341.   gimp_size_entry_set_unit (GIMP_SIZE_ENTRY (size_entry), GIMP_UNIT_PIXEL);
  2342.  
  2343.   /*  set the resolution to the image resolution  */
  2344.   gimp_size_entry_set_resolution (GIMP_SIZE_ENTRY (size_entry),
  2345.                                   0, xres, TRUE);
  2346.   gimp_size_entry_set_resolution (GIMP_SIZE_ENTRY (size_entry),
  2347.                                   1, yres, TRUE);
  2348.  
  2349.   /*  set the size (in pixels) that will be treated as 0% and 100%  */
  2350.   gimp_size_entry_set_size (GIMP_SIZE_ENTRY (size_entry),
  2351.                             0, 0.0, n);
  2352.   gimp_size_entry_set_size (GIMP_SIZE_ENTRY (size_entry),
  2353.                             1, 0.0, m);
  2354.  
  2355.   /*  set upper and lower limits (in pixels)  */
  2356.   gimp_size_entry_set_refval_boundaries (GIMP_SIZE_ENTRY (size_entry),
  2357.                                          0, n, MAX_WIDTH);
  2358.   gimp_size_entry_set_refval_boundaries (GIMP_SIZE_ENTRY (size_entry),
  2359.                                          1, m, MAX_HEIGHT);
  2360.  
  2361.   gtk_table_set_row_spacing (GTK_TABLE (size_entry), 0, 1);
  2362.   gtk_table_set_col_spacings (GTK_TABLE (size_entry), 6);
  2363.   gtk_table_set_col_spacing (GTK_TABLE (size_entry), 2, 12);
  2364.  
  2365.   /*  initialize the values  */
  2366.   gimp_size_entry_set_refval (GIMP_SIZE_ENTRY (size_entry), 0, n);
  2367.   gimp_size_entry_set_refval (GIMP_SIZE_ENTRY (size_entry), 1, m);
  2368.  
  2369.   /*  attach labels  */
  2370.   gimp_size_entry_attach_label (GIMP_SIZE_ENTRY (size_entry), "Width",
  2371.                                 0, 1, 0.0f);
  2372.   gimp_size_entry_attach_label (GIMP_SIZE_ENTRY (size_entry), "Height",
  2373.                                 0, 2, 0.0f);
  2374.  
  2375.   gtk_widget_show (size_entry);
  2376.   gtk_box_pack_start (GTK_BOX (main_hbox), size_entry, FALSE, FALSE, 0);
  2377.  
  2378.   /*  put a chain_button beside the size_entries  */
  2379.   aspect_chain = gimp_chain_button_new (GIMP_CHAIN_BOTTOM);
  2380.   gimp_chain_button_set_active (GIMP_CHAIN_BUTTON (aspect_chain), TRUE);
  2381.   gtk_table_attach_defaults (GTK_TABLE (size_entry),
  2382.                              aspect_chain, 1, 3, 2, 3);
  2383.   gtk_widget_show (aspect_chain);
  2384.  
  2385.   frame_label = gtk_label_new ("<b>Image Size</b>");
  2386.   gtk_widget_show (frame_label);
  2387.   gtk_frame_set_label_widget (GTK_FRAME (frame), frame_label);
  2388.   gtk_label_set_use_markup (GTK_LABEL (frame_label), TRUE);
  2389.  
  2390.   g_signal_connect (size_entry, "value-changed",
  2391.                     G_CALLBACK (scale_callback),
  2392.                     aspect_chain);
  2393.  
  2394.   gtk_widget_show (dialog);
  2395.   not_done = TRUE;
  2396.  
  2397.   do
  2398.   {
  2399.     run = gimp_dialog_run (GIMP_DIALOG (dialog));
  2400.     switch (run)
  2401.     {
  2402.       case RESPONSE_RESET:
  2403.       gimp_size_entry_set_refval (GIMP_SIZE_ENTRY (size_entry), 0, n);
  2404.       gimp_size_entry_set_refval (GIMP_SIZE_ENTRY (size_entry), 1, m);
  2405.       gimp_chain_button_set_active (GIMP_CHAIN_BUTTON (aspect_chain), TRUE);
  2406.       break;
  2407.  
  2408.       default:
  2409.       /* grab width and height from size entry boxes */
  2410.       svals.width  = (int) ((gimp_size_entry_get_refval
  2411.                            (GIMP_SIZE_ENTRY (size_entry), 0)) + 0.5);
  2412.       svals.height = (int) ((gimp_size_entry_get_refval
  2413.                            (GIMP_SIZE_ENTRY (size_entry), 1)) + 0.5);
  2414.       not_done = FALSE;
  2415.       break;
  2416.     }
  2417.   } while (not_done);
  2418.  
  2419.   gtk_widget_destroy (dialog);
  2420.  
  2421.   return run == GTK_RESPONSE_OK;
  2422. }
  2423.  
  2424. static void
  2425. scale_callback (GtkWidget *widget,
  2426.                 gpointer   data)
  2427. {
  2428.   gdouble new_width;
  2429.   gdouble new_height;
  2430.  
  2431.   new_width  = gimp_size_entry_get_refval (GIMP_SIZE_ENTRY (widget), 0);
  2432.   new_height = gimp_size_entry_get_refval (GIMP_SIZE_ENTRY (widget), 1);
  2433.  
  2434.   if (gimp_chain_button_get_active (GIMP_CHAIN_BUTTON (data)))
  2435.   {
  2436.     if (new_height != _ibfnbqh_height)
  2437.     {
  2438.       _ibfnbqh_height = new_height;
  2439.       _ibfnbqh_width = (int) (new_height / _ibfnbqh_m_over_n + 0.5);
  2440.       gimp_size_entry_set_refval (GIMP_SIZE_ENTRY (widget), 0, _ibfnbqh_width);
  2441.     }
  2442.     else
  2443.     {
  2444.       _ibfnbqh_width = new_width;
  2445.       _ibfnbqh_height = (int) (new_width * _ibfnbqh_m_over_n + 0.5);
  2446.       gimp_size_entry_set_refval (GIMP_SIZE_ENTRY (widget), 1, _ibfnbqh_height);
  2447.     }
  2448.   }
  2449.   else
  2450.   {
  2451.     _ibfnbqh_width  = new_width;
  2452.     _ibfnbqh_height = new_height;
  2453.   }
  2454. }

Paste is for source code and general debugging text.

Login or Register to edit, delete and keep track of your pastes and more.

Raw Paste

Login or Register to edit or fork this paste. It's free.