home *** CD-ROM | disk | FTP | other *** search
/ PC Pro 2002 April / pcpro0402.iso / essentials / graphics / Gimp / gimp-src-20001226.exe / src / gimp / unofficial-plug-ins / logconv / logconv.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-06-30  |  30.2 KB  |  1,065 lines
  1. /*
  2.  * logconv.c 1.2
  3.  *
  4.  * Copyright (C) 1998 Alessandro Baldoni
  5.  *
  6.  * This program is free software; you can redistribute it and/or modify
  7.  * it under the terms of the GNU General Public License as published by
  8.  * the Free Software Foundation; either version 2 of the License, or
  9.  * (at your option) any later version.
  10.  *
  11.  * This program is distributed in the hope that it will be useful,
  12.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13.  * MERCHANTABILITY of FITNESS FOR A PARTICULAR PURPOSE. See the
  14.  * GNU General Public License for more details.
  15.  *
  16.  * You should have received a copy of the GNU General Public License
  17.  * along with this program; if not, write to the Free Software
  18.  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  19.  */
  20. /*
  21.  * This is version 1.2
  22.  * LogConv for RGB images works on their luminance.
  23.  */
  24. #include <gtk/gtk.h>
  25. #include <libgimp/gimp.h>
  26. #include <math.h>
  27.  
  28. #ifdef __MSVCRT__
  29. #define isnan(x) _isnan(x)
  30. #endif
  31.  
  32. #define KD(sigma, ag, al) ((sigma * G_PI)/ (sqrt (pow (ag, 2.0) + pow(al, 2.0))))
  33. #define KS(sigma, kd, al) ((sigma * G_PI) / (al * kd))
  34. #define SG(sigma, ks) (sigma * sqrt (1.0 - (1.0 / pow (ks, 2.0))))
  35. #define SL(sigma, ks, kd) (sigma / (ks * kd))
  36.  
  37. typedef struct _entry
  38. {
  39.   gdouble Ag, Al;
  40. } entry;
  41.  
  42. static void ZoomImage (gdouble *, gdouble *, gint, gint, gint);
  43. static gdouble *LoGConvolution (GDrawable *, gdouble, gint);
  44. static gint LoGZero (GDrawable *, gdouble, gint, gint);
  45. static gdouble *MakeFilter (gdouble, gint, gint);
  46. static gint oddsupport (gdouble);
  47. static void query ();
  48. static void run (char *, int, GParam *, int *, GParam **);
  49. static gint LoGdialog (void);
  50. static void close_callback (GtkWidget *, gpointer);
  51. static void ok_callback (GtkWidget *, gpointer);
  52. static inline gdouble SobelX (gdouble *, gint, gint, gint);
  53. static inline gdouble SobelY (gdouble *, gint, gint, gint);
  54. static inline gdouble RobertsX (gdouble *, gint, gint, gint);
  55. static inline gdouble RobertsY (gdouble *, gint, gint, gint);
  56. static inline gdouble H1 (gdouble, gdouble);
  57. static inline gdouble H2 (gdouble, gdouble);
  58. static inline gdouble grad (gdouble, gdouble);
  59.  
  60. /*
  61.  * Precomputed values of Al e Ag for a given pa
  62.  */
  63. static entry Table1[11] =
  64. {
  65.   {3.554140, 4.205983},        /* pa =  0.0001 */
  66.   {3.402057, 4.061595},        /* pa =  0.0003 */
  67.   {3.227792, 3.896178},        /* pa =  0.0010 */
  68.   {3.060844, 3.737694},        /* pa =  0.0030 */
  69.   {2.867757, 3.554300},        /* pa =  0.0100 */
  70.   {2.712526, 3.406712},        /* pa =  0.0250 */
  71.   {2.461219, 3.167270},        /* pa =  0.1000 */
  72.   {2.244686, 2.960158},        /* pa =  0.3000 */
  73.   {1.984301, 2.709568},        /* pa =  1.0000 */
  74.   {1.718008, 2.450782},        /* pa =  3.0000 */
  75.   {1.378026, 2.115151}        /* pa = 10.0000 */
  76. };
  77.  
  78. typedef struct _LoGDialog
  79.   {
  80.     gint run, pa, pc1, pc2, which;
  81.     gdouble sd;
  82.   }
  83. LoGDialog;
  84.  
  85. static LoGDialog ldial = { FALSE, 0, 25, 60, 0, 2.0 };
  86.  
  87. static GtkWidget *toggle[11], *logentry, *ttoggle[3], *pc1entry, *pc2entry;
  88.  
  89. GPlugInInfo PLUG_IN_INFO =
  90. {
  91.   NULL,                /* init_proc */
  92.   NULL,                /* quit_proc */
  93.   query,            /* query_proc */
  94.   run,                /* run_proc */
  95. };
  96.  
  97. MAIN ()
  98.  
  99.      static void
  100.        query ()
  101. {
  102.   static GParamDef log_args[] =
  103.   {
  104.     {PARAM_INT32, "run_mode", "Interactive, non-interactive"},
  105.     {PARAM_IMAGE, "image", "Input image"},
  106.     {PARAM_DRAWABLE, "drawable", "Input drawable"},
  107.     {PARAM_INT32, "pa", "Selected pa (0..10)"},
  108.     {PARAM_FLOAT, "sigma", "Standard deviation"},
  109.     {PARAM_INT32, "type", "0: Standard LoG, 1: LoG with Roberts, 2: LoG with Sobel"}
  110.   };
  111.   static int nlog_args = sizeof (log_args) / sizeof (log_args[0]);
  112.  
  113.   gimp_install_procedure ("plug_in_LoG",
  114.               "Apply the LoG filter",
  115.               "",
  116.               "Alessandro Baldoni",
  117.               "Alessandro Baldoni",
  118.               "1998",
  119.               "<Image>/Filters/Edge-Detect/LoG",
  120.               "RGB,GRAY",
  121.               PROC_PLUG_IN,
  122.               nlog_args, 0,
  123.               log_args, NULL);
  124. }
  125.  
  126. static void
  127. run (char *name, int nparams, GParam * param, int *nreturn_vals,
  128.      GParam ** return_vals)
  129. {
  130.   static GParam values[1];
  131.   GRunModeType run_mode;
  132.   GDrawable *draw;
  133.  
  134.   run_mode = param[0].data.d_int32;
  135.  
  136.   *nreturn_vals = 1;
  137.   *return_vals = values;
  138.  
  139.   values[0].type = PARAM_STATUS;
  140.   values[0].data.d_status = STATUS_CALLING_ERROR;
  141.  
  142.   draw = gimp_drawable_get (param[2].data.d_int32);
  143.   switch (run_mode)
  144.     {
  145.     case RUN_INTERACTIVE:
  146.       gimp_get_data ("plug_in_LoG", &ldial);
  147.       if (LoGdialog ())
  148.     {
  149.       if (LoGZero (draw, ldial.sd, ldial.pa, ldial.which) != -1)
  150.         {
  151.           values[0].data.d_status = STATUS_SUCCESS;
  152.           gimp_displays_flush ();
  153.         }
  154.       else
  155.         values[0].data.d_status = STATUS_EXECUTION_ERROR;
  156.     }
  157.       break;
  158.     case RUN_NONINTERACTIVE:
  159.       if (LoGZero (draw, param[4].data.d_float, param[3].data.d_int32,
  160.            param[5].data.d_int32) != -1)
  161.     values[0].data.d_status = STATUS_SUCCESS;
  162.       else
  163.     values[0].data.d_status = STATUS_EXECUTION_ERROR;
  164.       break;
  165.     case RUN_WITH_LAST_VALS:
  166.       gimp_get_data ("plug_in_LoG", &ldial);
  167.       if (LoGZero (draw, ldial.sd, ldial.pa, ldial.which) != -1)
  168.     values[0].data.d_status = STATUS_SUCCESS;
  169.       else
  170.     values[0].data.d_status = STATUS_EXECUTION_ERROR;
  171.       break;
  172.     default:
  173.       break;
  174.     }
  175.   gimp_drawable_detach (draw);
  176. }
  177.  
  178. static gint
  179. LoGdialog ()
  180. {
  181.   GtkWidget *dlg, *button, *frame, *toggle_vbox, *hbox, *label, *hhbox;
  182.   GSList *group;
  183.   gchar **argv, *buf;
  184.   gint argc;
  185.   GtkTooltips *tips;
  186.  
  187.   argc = 1;
  188.   argv = g_new (gchar *, 1);
  189.   argv[0] = g_strdup ("save");
  190.   ldial.run = FALSE;
  191.  
  192.   gtk_init (&argc, &argv);
  193.   gtk_rc_parse (gimp_gtkrc ());
  194.  
  195.   tips = gtk_tooltips_new ();
  196.   dlg = gtk_dialog_new ();
  197.   gtk_window_set_title (GTK_WINDOW (dlg), "LoG filter");
  198.   gtk_window_position (GTK_WINDOW (dlg), GTK_WIN_POS_MOUSE);
  199.   gtk_window_set_wmclass (GTK_WINDOW (dlg), "LoG", "log");
  200.   gtk_signal_connect (GTK_OBJECT (dlg), "destroy",
  201.               (GtkSignalFunc) close_callback, NULL);
  202.   gtk_container_border_width (GTK_CONTAINER (dlg), 5);
  203.  
  204.   /*  Action area  */
  205.   button = gtk_button_new_with_label ("OK");
  206.   GTK_WIDGET_SET_FLAGS (button, GTK_CAN_DEFAULT);
  207.   gtk_signal_connect (GTK_OBJECT (button), "clicked",
  208.               (GtkSignalFunc) ok_callback, dlg);
  209.   gtk_box_pack_start (GTK_BOX (GTK_DIALOG (dlg)->action_area), button,
  210.               TRUE, TRUE, 0);
  211.   gtk_widget_grab_default (button);
  212.   gtk_widget_show (button);
  213.  
  214.   button = gtk_button_new_with_label ("Cancel");
  215.   GTK_WIDGET_SET_FLAGS (button, GTK_CAN_DEFAULT);
  216.   gtk_signal_connect_object (GTK_OBJECT (button), "clicked",
  217.                  (GtkSignalFunc) gtk_widget_destroy,
  218.                  GTK_OBJECT (dlg));
  219.   gtk_box_pack_start (GTK_BOX (GTK_DIALOG (dlg)->action_area), button,
  220.               TRUE, TRUE, 0);
  221.   gtk_widget_show (button);
  222.  
  223.   hbox = gtk_hbox_new (FALSE, 0);
  224.   gtk_box_pack_start (GTK_BOX (GTK_DIALOG (dlg)->vbox), hbox, FALSE, TRUE, 0);
  225.   gtk_widget_show (hbox);
  226.  
  227.   /*  compression  */
  228.   frame = gtk_frame_new ("Allowable PA");
  229.   gtk_frame_set_shadow_type (GTK_FRAME (frame), GTK_SHADOW_ETCHED_IN);
  230.   gtk_container_border_width (GTK_CONTAINER (frame), 10);
  231.   gtk_box_pack_start (GTK_BOX (hbox), frame, FALSE, TRUE, 0);
  232.   toggle_vbox = gtk_table_new (6, 2, FALSE);
  233.   gtk_container_border_width (GTK_CONTAINER (toggle_vbox), 5);
  234.   gtk_container_add (GTK_CONTAINER (frame), toggle_vbox);
  235.  
  236.   group = NULL;
  237.   toggle[0] = gtk_radio_button_new_with_label (group, " 0.0001%");
  238.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[0]));
  239.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[0], 0, 1, 0, 1);
  240.   gtk_widget_show (toggle[0]);
  241.  
  242.   toggle[1] = gtk_radio_button_new_with_label (group, " 0.0003%");
  243.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[1]));
  244.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[1], 0, 1, 1, 2);
  245.   gtk_widget_show (toggle[1]);
  246.  
  247.   toggle[2] = gtk_radio_button_new_with_label (group, " 0.0010%");
  248.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[2]));
  249.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[2], 0, 1, 2, 3);
  250.   gtk_widget_show (toggle[2]);
  251.  
  252.   toggle[3] = gtk_radio_button_new_with_label (group, " 0.0030%");
  253.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[3]));
  254.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[3], 0, 1, 3, 4);
  255.   gtk_widget_show (toggle[3]);
  256.  
  257.   toggle[4] = gtk_radio_button_new_with_label (group, " 0.0100%");
  258.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[4]));
  259.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[4], 0, 1, 4, 5);
  260.   gtk_widget_show (toggle[4]);
  261.  
  262.   toggle[5] = gtk_radio_button_new_with_label (group, " 0.0250%");
  263.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[5]));
  264.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[5], 0, 1, 5, 6);
  265.   gtk_widget_show (toggle[5]);
  266.  
  267.   toggle[6] = gtk_radio_button_new_with_label (group, " 0.1000%");
  268.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[6]));
  269.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[6], 1, 2, 0, 1);
  270.   gtk_widget_show (toggle[6]);
  271.  
  272.   toggle[7] = gtk_radio_button_new_with_label (group, " 0.3000%");
  273.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[7]));
  274.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[7], 1, 2, 1, 2);
  275.   gtk_widget_show (toggle[7]);
  276.  
  277.   toggle[8] = gtk_radio_button_new_with_label (group, " 1.0000%");
  278.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[8]));
  279.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[8], 1, 2, 2, 3);
  280.   gtk_widget_show (toggle[8]);
  281.  
  282.   toggle[9] = gtk_radio_button_new_with_label (group, " 3.0000%");
  283.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[9]));
  284.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[9], 1, 2, 3, 4);
  285.   gtk_widget_show (toggle[9]);
  286.  
  287.   toggle[10] = gtk_radio_button_new_with_label (group, "10.0000%");
  288.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (toggle[10]));
  289.   gtk_table_attach_defaults (GTK_TABLE (toggle_vbox), toggle[10], 1, 2, 4, 5);
  290.   gtk_widget_show (toggle[10]);
  291.  
  292.   gtk_toggle_button_set_state (GTK_TOGGLE_BUTTON (toggle[ldial.pa]),
  293.                    TRUE);
  294.  
  295.   gtk_widget_show (toggle_vbox);
  296.   gtk_widget_show (frame);
  297.  
  298.   /* LoG type */
  299.   frame = gtk_frame_new ("LoG type");
  300.   gtk_frame_set_shadow_type (GTK_FRAME (frame), GTK_SHADOW_ETCHED_IN);
  301.   gtk_container_border_width (GTK_CONTAINER (frame), 10);
  302.   gtk_box_pack_start (GTK_BOX (hbox), frame, FALSE, TRUE, 0);
  303.   toggle_vbox = gtk_vbox_new (FALSE, 5);
  304.   gtk_container_border_width (GTK_CONTAINER (toggle_vbox), 5);
  305.   gtk_container_add (GTK_CONTAINER (frame), toggle_vbox);
  306.  
  307.   group = NULL;
  308.   ttoggle[0] = gtk_radio_button_new_with_label (group, "Standard LoG");
  309.   gtk_tooltips_set_tip (tips, ttoggle[0], "Do not apply any gradient or thresholding", "");
  310.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (ttoggle[0]));
  311.   gtk_box_pack_start (GTK_BOX (toggle_vbox), ttoggle[0], FALSE, FALSE, 0);
  312.   gtk_widget_show (ttoggle[0]);
  313.  
  314.   ttoggle[1] = gtk_radio_button_new_with_label (group, "LoG with Roberts");
  315.   gtk_tooltips_set_tip (tips, ttoggle[1], "Apply a Roberts gradient and a threshold", "");
  316.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (ttoggle[1]));
  317.   gtk_box_pack_start (GTK_BOX (toggle_vbox), ttoggle[1], FALSE, FALSE, 0);
  318.   gtk_widget_show (ttoggle[1]);
  319.  
  320.   ttoggle[2] = gtk_radio_button_new_with_label (group, "LoG with Sobel");
  321.   gtk_tooltips_set_tip (tips, ttoggle[2], "Apply a Sobel gradient and a threshold", "");
  322.   group = gtk_radio_button_group (GTK_RADIO_BUTTON (ttoggle[2]));
  323.   gtk_box_pack_start (GTK_BOX (toggle_vbox), ttoggle[2], FALSE, FALSE, 0);
  324.   gtk_widget_show (ttoggle[2]);
  325.  
  326.   gtk_toggle_button_set_state (GTK_TOGGLE_BUTTON (ttoggle[ldial.which]),
  327.                    TRUE);
  328.   gtk_widget_show (toggle_vbox);
  329.   gtk_widget_show (frame);
  330.  
  331.   /* Standard deviation */
  332.   hhbox = gtk_hbox_new (FALSE, 0);
  333.   gtk_box_pack_start (GTK_BOX (GTK_DIALOG (dlg)->vbox), hhbox, FALSE, TRUE, 0);
  334.   gtk_widget_show (hhbox);
  335.  
  336.   label = gtk_label_new ("Standard deviation: ");
  337.   gtk_box_pack_start (GTK_BOX (hhbox), label, FALSE, TRUE, 0);
  338.   gtk_widget_show (label);
  339.  
  340.   buf = g_new (gchar, 10);
  341.   sprintf (buf, "%f", ldial.sd);
  342.   logentry = gtk_entry_new ();
  343.   gtk_tooltips_set_tip (tips, logentry, "This is the standard deviation: the smaller it is, the finer details you'll get. The bigger it is, the coarser details you'll get. If this is smaller than 2.0, you must choose a big (3.0, 10.0) PA", "");
  344.   gtk_box_pack_start (GTK_BOX (hhbox), logentry, FALSE, TRUE, 0);
  345.   gtk_entry_set_text (GTK_ENTRY (logentry), buf);
  346.   g_free (buf);
  347.   gtk_widget_show (logentry);
  348.  
  349.   /* pc1 */
  350.   label = gtk_label_new ("Thresholds used after gradient computation ");
  351.   gtk_misc_set_alignment (GTK_MISC (label), 0.0, 0.5);
  352.   gtk_box_pack_start (GTK_BOX (GTK_DIALOG (dlg)->vbox), label, FALSE, TRUE, 0);
  353.   gtk_widget_show (label);
  354.  
  355.   hhbox = gtk_hbox_new (FALSE, 0);
  356.   gtk_box_pack_start (GTK_BOX (GTK_DIALOG (dlg)->vbox), hhbox, FALSE, TRUE, 0);
  357.   gtk_widget_show (hhbox);
  358.  
  359.   label = gtk_label_new ("PC1: ");
  360.   gtk_box_pack_start (GTK_BOX (hhbox), label, FALSE, TRUE, 0);
  361.   gtk_widget_show (label);
  362.  
  363.   buf = g_new (gchar, 10);
  364.   sprintf (buf, "%d", ldial.pc1);
  365.   pc1entry = gtk_entry_new ();
  366.   gtk_tooltips_set_tip (tips, pc1entry, "If you choosed to apply a gradient, values below this threshold will be skipped", "");
  367.   gtk_box_pack_start (GTK_BOX (hhbox), pc1entry, FALSE, TRUE, 0);
  368.   gtk_entry_set_text (GTK_ENTRY (pc1entry), buf);
  369.   g_free (buf);
  370.   gtk_widget_show (pc1entry);
  371.  
  372.   /* pc2 */
  373.   hhbox = gtk_hbox_new (FALSE, 0);
  374.   gtk_box_pack_start (GTK_BOX (GTK_DIALOG (dlg)->vbox), hhbox, FALSE, TRUE, 0);
  375.   gtk_widget_show (hhbox);
  376.  
  377.   label = gtk_label_new ("PC2: ");
  378.   gtk_box_pack_start (GTK_BOX (hhbox), label, FALSE, TRUE, 0);
  379.   gtk_widget_show (label);
  380.  
  381.   buf = g_new (gchar, 10);
  382.   sprintf (buf, "%d", ldial.pc2);
  383.   pc2entry = gtk_entry_new ();
  384.   gtk_tooltips_set_tip (tips, pc2entry, "If you choosed to apply a gradient, values exceeding this threshold will be skipped", "");
  385.   gtk_box_pack_start (GTK_BOX (hhbox), pc2entry, FALSE, TRUE, 0);
  386.   gtk_entry_set_text (GTK_ENTRY (pc2entry), buf);
  387.   g_free (buf);
  388.   gtk_widget_show (pc2entry);
  389.  
  390.   gtk_widget_show (dlg);
  391.  
  392.   gtk_main ();
  393.   gdk_flush ();
  394.  
  395.   return ldial.run;
  396. }
  397.  
  398. static void
  399. close_callback (GtkWidget * widget, gpointer data)
  400. {
  401.   gtk_main_quit ();
  402. }
  403.  
  404. static void
  405. ok_callback (GtkWidget * widget, gpointer data)
  406. {
  407.   gint i;
  408.  
  409.   for (i = 0; i < 11; i++)
  410.     if (GTK_TOGGLE_BUTTON (toggle[i])->active)
  411.       ldial.pa = i;
  412.   for (i = 0; i < 3; i++)
  413.     if (GTK_TOGGLE_BUTTON (ttoggle[i])->active)
  414.       ldial.which = i;
  415.   ldial.sd = atof (gtk_entry_get_text (GTK_ENTRY (logentry)));
  416.   ldial.pc1 = atoi (gtk_entry_get_text (GTK_ENTRY (pc1entry)));
  417.   ldial.pc2 = atoi (gtk_entry_get_text (GTK_ENTRY (pc2entry)));
  418.   ldial.run = TRUE;
  419.   gimp_set_data ("plug_in_LoG", &ldial, sizeof (LoGDialog));
  420.   gtk_widget_destroy (GTK_WIDGET (data));
  421. }
  422.  
  423. /*
  424.  * Computes the convolution of draw with a LoG of constant sigma and 
  425.  * aliasing pa.
  426.  */
  427. static gdouble *
  428. LoGConvolution (GDrawable * draw, gdouble sigma, gint pa)
  429. {
  430.   gint i, j, k, addr, MSL, W, H, W1, H1, s2, s3, kd, MG, i1, j1, j2, gap;
  431.   gdouble *temp, *temp2, *h1, *h2, *g, h1x, h1y, h2x, h2y, *temp1, 
  432.     Al, Ag, ks, sigmal, sigmag, *temp3, t;
  433.   GPixelRgn src;
  434.   guchar *buf;
  435.  
  436.   W = gimp_drawable_width (draw->id);
  437.   H = gimp_drawable_height (draw->id);
  438.   gap = (gimp_drawable_type (draw->id) == RGB_IMAGE) ? 3 : 1;
  439.   buf = g_new (guchar, W * H * gap);
  440.   gimp_pixel_rgn_init (&src, draw, 0, 0, W, H, FALSE, FALSE);
  441.   gimp_pixel_rgn_get_rect (&src, buf, 0, 0, W, H);
  442.   gimp_progress_init ("Computing LoG...");
  443.   /* 
  444.    * 1. Operator design 
  445.    */
  446.   /* Cutoff constants */
  447.   Ag = Table1[pa].Ag;
  448.   Al = Table1[pa].Al;
  449.   /* Decimation factor */
  450.   kd = (gint) floor (KD (sigma, Ag, Al));
  451.   /* Reconstruction constant */
  452.   ks = KS (sigma, (gdouble) kd, Al);
  453.   /* Gaussian space constant */
  454.   sigmag = SG (sigma, ks);
  455.   /* LoG space constant */
  456.   sigmal = SL (sigma, ks, kd);
  457.   if ((isnan (sigmag) != 0) || (isnan (sigmal) != 0))
  458.     return (gdouble *) NULL;
  459.   /* LoG MSL width */
  460.   MSL = oddsupport (sigmal);
  461.   /* LoG filters computation */
  462.   h1 = MakeFilter (sigmal, MSL, 1);
  463.   h2 = MakeFilter (sigmal, MSL, 2);
  464.   /* Gaussian MSL width */
  465.   MG = 6 * (gint) sigmag;
  466.   if (MG % 2 == 0)
  467.     MG++;
  468.   /* Gaussian filter computation */
  469.   g = MakeFilter (sigmag, MG, 3);
  470.   /*
  471.    * 2. Convolution process
  472.    */
  473.   /* Gaussian convolution and decimation */
  474.   W1 = (gint) ceil ((gdouble) W / (gdouble) kd);
  475.   H1 = (gint) ceil ((gdouble) H / (gdouble) kd);
  476.   s3 = MG / 2;
  477.   temp1 = g_new (gdouble, W1 * (H + s3 * 2));
  478.   temp = g_new (gdouble, (W + s3 * 2) * H);
  479.   i1 = W + s3 * 2;
  480.   if (gap == 1)                 /* Gray images */
  481.     for (i = 0; i < H; i++)
  482.       {
  483.     j1 = i * i1;
  484.     for (j = 0; j < W; j++)
  485.       temp[(j + s3) + j1] = ((gdouble) buf[j + i * W]);
  486.       }
  487.   else                          /* RGB images */
  488.     for (i = 0; i < H; i++)
  489.       {
  490.     j1 = i * i1;
  491.     for (j2 = 0, j = 0; j < W; j++, j2 += 3)
  492.       /* Luminance */
  493.       temp[(j + s3) + j1] = (gdouble) (buf[j2 + i * W * 3] * 0.299 + 
  494.                        buf[j2 + 1 + i * W * 3] * 0.587 +
  495.                        buf[j2 + 2 + i * W * 3] * 0.114);
  496.       }
  497.   g_free (buf);
  498.   gimp_progress_update (1.0 / 7.0);
  499.   /* DC-padding */
  500.   for (i = 0; i < H; i++)
  501.     {                /* Left */
  502.       j1 = i * (W + s3 * 2);
  503.       i1 = s3 + j1;
  504.       for (j = 0; j < s3; j++)
  505.     temp[j + j1] = temp[i1];
  506.     }
  507.   for (i = 0; i < H; i++)
  508.     {                /* Right */
  509.       j1 = i * (W + s3 * 2);
  510.       i1 = (W + s3 - 1) + j1;
  511.       for (j = W + s3; j < (W + s3 * 2); j++)
  512.     temp[j + j1] = temp[i1];
  513.     }
  514.   for (i = 0, i1 = s3; i < H; i++, i1++)    /* Row convolution */
  515.     for (j = s3, j1 = 0; j < W + s3; j += kd, j1++)
  516.       {
  517.     t = 0.0;
  518.     for (k = 0; k < MG; k++)
  519.       t += temp[(j - s3 + k) + i * (W + s3 * 2)] * g[k];
  520.     temp1[j1 + i1 * W1] = t;
  521.       }
  522.   g_free (temp);
  523.   gimp_progress_update (2.0 / 7.0);
  524.   /* DC-padding */
  525.   i1 = s3 * W1;
  526.   for (i = 0; i < s3; i++)
  527.     {                /* Up */
  528.       j1 = i * W1;
  529.       for (j = 0; j < W1; j++)
  530.     temp1[j + j1] = temp1[j + i1];
  531.     }
  532.   i1 = (H + s3 - 1) * W1;
  533.   for (i = H + s3; i < (H + s3 * 2); i++)
  534.     {                /* Down */
  535.       j1 = i * W1;
  536.       for (j = 0; j < W1; j++)
  537.     temp1[j + j1] = temp1[j + i1];
  538.     }
  539.   temp2 = g_new (gdouble, W1 * H1);
  540.   for (i = s3, i1 = 0; i < H + s3; i += kd, i1++)     /* Column convolution */
  541.     for (j = 0; j < W1; j++)
  542.       {
  543.     t = 0.0;
  544.     for (k = 0; k < MG; k++)
  545.       t += temp1[j + (i - s3 + k) * W1] * g[k];
  546.     temp2[j + i1 * W1] = t;
  547.       }
  548.   g_free (temp1);
  549.   gimp_progress_update (3.0 / 7.0);
  550.   /* LoG convolution */
  551.   s2 = MSL / 2;
  552.   temp = g_new (gdouble, (W1 + s2 * 2) * (H1 + s2 * 2));
  553.   i1 = W1 + s2 * 2;
  554.   for (i = 0; i < H1; i++)
  555.     for (j = 0; j < W1; j++)
  556.       temp[(j + s2) + (i + s2) * i1] = temp2[j + i * W1];
  557.   gimp_progress_update (4.0 / 7.0);
  558.   /* DC-padding */
  559.   i1 = W1 + s2 * 2;
  560.   for (i = s2; i < H1 + s2; i++)
  561.     {                /* Left */
  562.       j1 = i * i1;
  563.       for (j = 0; j < s2; j++)
  564.     temp[j + j1] = temp[s2 + j1];
  565.     }
  566.   for (i = s2; i < H1 + s2; i++)
  567.     {                /* Right */
  568.       j1 = i * i1;
  569.       for (j = W1 + s2; j < (W1 + s2 * 2); j++)
  570.     temp[j + j1] = temp[(W1 + s2 - 1) + j1];
  571.     }
  572.   g_free (temp2);
  573.   temp1 = g_new (gdouble, W1 * (H1 + s2 * 2));
  574.   temp3 = g_new (gdouble, W1 * (H1 + s2 * 2));
  575.   for (i = s2, i1 = 0; i < H1 + s2; i++, i1++)
  576.     {                /* Row convolution */
  577.       for (j = s2, j1 = 0; j < W1 + s2; j++, j1++)
  578.     {
  579.       addr = i1 * W1;
  580.       h1x = h2x = 0.0;
  581.       for (k = 0; k < MSL; k++)
  582.         {
  583.           h1x += temp[(j - s2 + k) + i * (W1 + s2 * 2)] * h1[k];
  584.           h2x += temp[(j - s2 + k) + i * (W1 + s2 * 2)] * h2[k];
  585.         }
  586.       temp1[j1 + i * W1] = h1x;
  587.       temp3[j1 + i * W1] = h2x;
  588.     }
  589.     }
  590.   g_free (temp);
  591.   gimp_progress_update (5.0 / 7.0);
  592.   /* DC-padding */
  593.   i1 = s2 * W1;
  594.   for (i = 0; i < s2; i++)
  595.     {                /* Up */
  596.       j1 = i * W1;
  597.       for (j = 0; j < W1; j++)
  598.     {
  599.       temp1[j + j1] = temp1[j + i1];
  600.       temp3[j + j1] = temp3[j + i1];
  601.     }
  602.     }
  603.   i1 = (H1 + s2 - 1) * W1;
  604.   for (i = H1 + s2; i < (H1 + s2 * 2); i++)
  605.     {                /* Down */
  606.       j1 = i * W1;
  607.       for (j = 0; j < W1; j++)
  608.     {
  609.       temp1[j + j1] = temp1[j + i1];
  610.       temp3[j + j1] = temp3[j + i1];
  611.     }
  612.     }
  613.   temp = g_new (gdouble, W1 * H1);
  614.   for (i = s2, i1 = 0; i < H1 + s2; i++, i1++)
  615.     {                /* Column convolution */
  616.       for (j = 0; j < W1; j++)
  617.     {
  618.       h1y = h2y = 0.0;
  619.       for (k = 0; k < MSL; k++)
  620.         {
  621.           h1y += temp3[j + (i - s2 + k) * W1] * h1[k];
  622.           h2y += temp1[j + (i - s2 + k) * W1] * h2[k];
  623.         }
  624.       temp[j + i1 * W1] = h1y + h2y;
  625.     }
  626.     }
  627.   g_free (temp1);
  628.   g_free (temp3);
  629.   gimp_progress_update (6.0 / 7.0);
  630.   /* Expansion */
  631.   temp1 = g_new (gdouble, W * H);
  632.   if (kd > 1)
  633.     ZoomImage (temp, temp1, kd, W, H);
  634.   else
  635.     for (i = 0; i < H * W; i++)
  636.       temp1[i] = temp[i];
  637.   gimp_progress_update (7.0 / 7.0);
  638.   g_free (temp);
  639.   g_free (h2);
  640.   g_free (h1);
  641.   g_free (g);
  642.   return temp1;
  643. }
  644.  
  645. /*
  646.  * Computes zero crossings for an image convoluted with a LoG of constant 
  647.  * sigma and aliasing pa.
  648.  * If gradient is specified, it performs a gradient thresholding as in
  649.  * A. Basu et al. (1995).
  650.  * gradient = 1 -> Roberts, = 2 -> Sobel
  651.  */
  652. static gint
  653. LoGZero (GDrawable * draw, gdouble sigma, gint pa, gint gradient)
  654. {
  655.   gdouble *temp1, pix0, pix1, pix2, pix3, pix4, gradx, grady, *gradval, 
  656.     max, pc1, pc2;
  657.   gint W, H, addr, i, j, gap, j1, i1;
  658.   gdouble (*funcx) (gdouble *, gint, gint, gint),
  659.     (*funcy) (gdouble *, gint, gint, gint);
  660.   GPixelRgn src;
  661.   guchar *buf;
  662.   gboolean do_zero;
  663.  
  664.   gradval = (gdouble *) NULL;
  665.   max = 0.0;
  666.   funcx = NULL;
  667.   funcy = NULL;
  668.   W = gimp_drawable_width (draw->id);
  669.   H = gimp_drawable_height (draw->id);
  670.   gap = (gimp_drawable_type (draw->id) == RGB_IMAGE) ? 3 : 1;
  671.   if ((temp1 = LoGConvolution (draw, sigma, pa)) == (gdouble *) NULL)
  672.     {
  673.       g_message ("LoG: Please choose a bigger PA.");
  674.       return -1;
  675.     }
  676.   /*
  677.    * Zero crossing detection and gradient computation.
  678.    */
  679.   if (gradient != 0)
  680.     {
  681.       max = 0.0;
  682.       gradval = g_new (gdouble, W * H);
  683.       memset (gradval, 0, sizeof (gdouble) * W * H);
  684.       if (gradient == 1)
  685.     {
  686.       funcx = RobertsX;
  687.       funcy = RobertsY;
  688.     }
  689.       else
  690.     {
  691.       funcx = SobelX;
  692.       funcy = SobelY;
  693.     }
  694.     }
  695.   buf = g_new (guchar, W * H * gap);
  696.   memset (buf, 255, sizeof (guchar) * W * H * gap);
  697.   gimp_progress_init ("Zero crossing...");
  698.   for (i = 1; i < H - 1; i++)
  699.     {
  700.       for (j1 = 3, j = 1; j < W - 1; j++, j1 += 3)
  701.     {
  702.       do_zero = FALSE;
  703.       addr = j + i * W;
  704.       /*
  705.        * Chessboard metric
  706.        *     1
  707.        *   2 0 3
  708.        *     4
  709.        */
  710.       pix0 = temp1[addr];
  711.       pix1 = temp1[j + (i - 1) * W];
  712.       pix2 = temp1[j - 1 + i * W];
  713.       pix3 = temp1[j + 1 + i * W];
  714.       pix4 = temp1[j + (i + 1) * W];
  715.       /*
  716.        * Zero crossing test by Simon A. J. Winder (c) 1994
  717.        */
  718.       if (pix0 > 0.0 && 
  719.           (pix1 < 0.0 || pix2 < 0.0 || pix3 < 0.0 || pix4 < 0.0))
  720.         do_zero = TRUE;
  721.       if (pix0 == 0.0)
  722.         {
  723.           if ((pix1 > 0.0 && pix4 < 0.0) || (pix1 < 0.0 && pix4 > 0.0) ||
  724.           (pix2 > 0.0 && pix3 < 0.0) || (pix2 < 0.0 && pix3 > 0.0))
  725.         do_zero = TRUE;
  726.           else
  727.         {
  728.           pix1 = temp1[j - 1 + (i + 1) * W];
  729.           pix2 = temp1[j + 1 + (i + 1) * W];
  730.           pix3 = temp1[j - 1 + (i - 1) * W];
  731.           pix4 = temp1[j + 1 + (i - 1) * W];
  732.           if ((pix1 > 0.0 && pix4 < 0.0) || 
  733.               (pix1 < 0.0 && pix4 > 0.0) ||
  734.               (pix2 > 0.0 && pix3 < 0.0) || (pix2 < 0.0 && pix3 > 0.0))
  735.             do_zero = TRUE;
  736.         }
  737.         }
  738.       if (do_zero)
  739.         {
  740.           if (gradient != 0)
  741.         {
  742.           gradx = funcx (temp1, j, i, W);
  743.           grady = funcy (temp1, j, i, W);
  744.           gradval[addr] = fabs (gradx) + fabs (grady);
  745.           if (max < gradval[addr])
  746.             max = gradval[addr];
  747.         }
  748.           else
  749.         {
  750.           if (gap == 3)
  751.             buf[j1 + i * W * 3] =  buf[j1 + 1 + i * W * 3] = 
  752.               buf[j1 + 2 + i * W * 3] = 0;
  753.           else
  754.             buf[addr] = 0;
  755.         }
  756.         }
  757.     }
  758.       gimp_progress_update ((gdouble) i / (gdouble) (H - 1));
  759.     }
  760.   if (gradient != 0)
  761.     {
  762.       /*
  763.        * Gradient weighting
  764.        */
  765.       pc1 = (max * (gdouble) ldial.pc1) / 100.0;    /* 25% */
  766.       pc2 = (max * (gdouble) ldial.pc2) / 100.0;    /* 60% */
  767.       /*
  768.        * We skip all pixels below the 25% gradient threshold
  769.        */
  770.       if (gap == 1)          /* Gray images */
  771.     {
  772.       for (i = 0; i < W * H; i++)
  773.         if ((gradval[i] > pc1) && (gradval[i] < pc2))
  774.           buf[i] = 0;
  775.     }
  776.       else                   /* RGB images */
  777.     for (i1 = 0, i = 0; i < W * H; i++, i1 += 3)
  778.       if ((gradval[i] > pc1) && (gradval[i] < pc2))
  779.         buf[i1] = buf[i1 + 1] = buf[i1 + 2] = 0;
  780.       g_free (gradval);
  781.     }
  782.   gimp_pixel_rgn_init (&src, draw, 0, 0, W, H, TRUE, TRUE);
  783.   gimp_pixel_rgn_set_rect (&src, buf, 0, 0, W, H);
  784.   gimp_drawable_flush (draw);
  785.   gimp_drawable_merge_shadow (draw->id, TRUE);
  786.   gimp_drawable_update (draw->id, 0, 0, W, H);
  787.   g_free (buf);
  788.   g_free (temp1);
  789.   return 0;
  790. }
  791.  
  792. /*
  793.  * Computes the filter's support. It must be odd.
  794.  */
  795. static gint
  796. oddsupport (gdouble sigma)
  797. {
  798.   gint support;
  799.  
  800.   support = 4 * 2 * G_SQRT2 * sigma;
  801.   if ((support % 2) == 0)
  802.     return (support + 1);
  803.   else
  804.     return support;
  805. }
  806.  
  807. /*
  808.  * Creates filters for LoG decomposition
  809.  */
  810. static gdouble *
  811. MakeFilter (gdouble sigma, gint support, gint which)
  812. {
  813.   gint x;
  814.   gdouble *filter;
  815.   gdouble (*func) (gdouble, gdouble);
  816.  
  817.   func = NULL;
  818.   switch (which)
  819.     {
  820.     case 1:
  821.       func = H1;
  822.       break;
  823.     case 2:
  824.       func = H2;
  825.       break;
  826.     case 3:
  827.       func = grad;
  828.       break;
  829.     }
  830.   filter = g_new (gdouble, support);
  831.   for (x = 0; x <= support / 2; x++)
  832.     filter[support / 2 + x] = filter[support / 2 - x] = func ((gdouble) x,
  833.                                   sigma);
  834.   return filter;
  835. }
  836.  
  837. static inline gdouble
  838. SobelX (gdouble * temp1, gint j, gint i, gint W)
  839. {
  840.   return (-temp1[(j - 1) + (i - 1) * W] - 2.0 * temp1[j + (i - 1) * W] -
  841.       temp1[(j + 1) + (i - 1) * W] + temp1[(j - 1) + (i + 1) * W] +
  842.       2.0 * temp1[j + (i + 1) * W] + temp1[(j + 1) + (i + 1) * W]);
  843. }
  844.  
  845. static inline gdouble
  846. SobelY (gdouble * temp1, gint j, gint i, gint W)
  847. {
  848.   return (-temp1[(j - 1) + (i - 1) * W] - 2.0 * temp1[(j - 1) + i * W] -
  849.       temp1[(j - 1) + (i + 1) * W] + temp1[(j + 1) + (i - 1) * W] +
  850.       2.0 * temp1[(j + 1) + i * W] + temp1[(j + 1) + (i + 1) * W]);
  851. }
  852.  
  853. static inline gdouble
  854. RobertsX (gdouble * temp1, gint j, gint i, gint W)
  855. {
  856.   return (temp1[j + i * W] - temp1[(j + 1) + (i + 1) * W]);
  857. }
  858.  
  859. static inline gdouble
  860. RobertsY (gdouble * temp1, gint j, gint i, gint W)
  861. {
  862.   return (temp1[(j + 1) + i * W] - temp1[j + (i + 1) * W]);
  863. }
  864.  
  865. static inline gdouble
  866. H1 (gdouble xi, gdouble sigma)
  867. {
  868.   return (1.0 / (sqrt (2.0 * G_PI) * pow (sigma, 2.0))) *
  869.     (1.0 - (pow (xi, 2.0) / pow (sigma, 2.0))) *
  870.     exp (-pow (xi, 2.0) / (2.0 * pow (sigma, 2.0)));
  871. }
  872.  
  873. static inline gdouble
  874. H2 (gdouble xi, gdouble sigma)
  875. {
  876.   return (1.0 / (sqrt (2.0 * G_PI) * pow (sigma, 2.0))) *
  877.     exp (-pow (xi, 2.0) / (2.0 * pow (sigma, 2.0)));
  878. }
  879.  
  880. static inline gdouble
  881. grad (gdouble x, gdouble sigma)
  882. {
  883.   return ((1.0 / (sqrt (2.0 * G_PI) * sigma)) *
  884.       exp (-pow (x, 2.0) / (2 * pow (sigma, 2.0))));
  885. }
  886.  
  887. /*
  888.  * A C-port from LAPACK...
  889.  */
  890. static inline void 
  891. dgefac (gdouble a[4][4], gdouble b[4])
  892. {
  893.   gdouble t, dmax, ipvt[4];
  894.   gint k, i, l, j;
  895.  
  896.   /*
  897.    * Gaussian elimination with partial pivoting 
  898.    */
  899.   for (k = 0; k < 3; k++)
  900.     {
  901.       /*
  902.        * Look for pivot l
  903.        */
  904.       l = k;
  905.       dmax = fabs (a[k][k]);
  906.       for (i = k + 1; i < 4; i++)
  907.     if (fabs (a[i][k]) > dmax)
  908.       {
  909.         dmax = fabs (a[i][k]);
  910.         l = i;
  911.       }
  912.       ipvt[k] = l;
  913.       /*
  914.        * Swap if needed
  915.        */
  916.       if (l != k)
  917.     {
  918.       t = a[l][k];
  919.       a[l][k] = a[k][k];
  920.       a[k][k] = t;
  921.     }
  922.       /*
  923.        * Multipliers computation
  924.        */
  925.       t = -1.0 / a[k][k];
  926.       for (i = k + 1; i < 4; i++)
  927.     a[i][k] *= t;
  928.       /*
  929.        * Row elimination
  930.        */
  931.       for (j = k + 1; j < 4; j++)
  932.     {
  933.       t = a[l][j];
  934.       if (l != k)
  935.         {
  936.           a[l][j] = a[k][j];
  937.           a[k][j] = t;
  938.         }
  939.       for (i = k + 1; i < 4; i++)
  940.         a[i][j] += (t * a[i][k]);
  941.     }
  942.     }
  943.   ipvt[3] = 3;
  944.   /*
  945.    * Solving L y = b
  946.    */
  947.   for (k = 0; k < 3; k++)
  948.     {
  949.       l = ipvt[k];
  950.       t = b[l];
  951.       if (l != k)
  952.     {
  953.       b[l] = b[k];
  954.       b[k] = t;
  955.     }
  956.       for (i = k + 1; i < 4; i++)
  957.     b[i] += (t * a[i][k]);
  958.     }
  959.   /*
  960.    * Then U x = y
  961.    */
  962.   for (k = 3; k >= 0; k--)
  963.     {
  964.       b[k] /= a[k][k];
  965.       t = -b[k];
  966.       for (i = 0; i < k; i++)
  967.     b[i] += (a[i][k] * t);
  968.     }
  969. }
  970.  
  971. /*
  972.  * Upsampling through bilinear interpolation.
  973.  */
  974. static void
  975. ZoomImage (gdouble * src, gdouble * dest, gint DF, gint width, gint height)
  976. {
  977.   gdouble b[4], *temp;
  978.   gdouble a[4][4];
  979.   gint i, j, l, k, W, H, addr, addr2;
  980.  
  981.   /*
  982.    * Buffer initialization
  983.    */
  984.   temp = g_new (gdouble, (width + DF) * (height + DF));
  985.   memset (temp, 0, sizeof (gdouble) * (width + DF) * (height + DF));
  986.   W = (gint) ceil ((gdouble) width / (gdouble) DF);
  987.   H = (gint) ceil ((gdouble) height / (gdouble) DF);
  988.   for (i = 0, k = 0; k < H; i += DF, k++)
  989.     for (j = 0, l = 0; l < W; j += DF, l++)
  990.       temp[j + i * (width + DF)] = src[l + k * W];
  991.   /*
  992.    * DC-padding
  993.    */
  994.   k = W * DF;
  995.   l = (W - 1) * DF;
  996.   for (i = 0; i < height + DF; i += DF)
  997.     temp[k + i * (width + DF)] = temp[l + i * (width + DF)];
  998.   k = H * DF * (width + DF);
  999.   l = (H - 1) * DF * (width + DF);
  1000.   for (j = 0; j < width + DF; j += DF)
  1001.     temp[j + k] = temp[j + l];
  1002.   W = width + DF;
  1003.   /*
  1004.    * Row bilinear interpolation 
  1005.    */
  1006.   for (i = 0; i < height; i += DF)
  1007.     {
  1008.       addr = i * W;
  1009.       addr2 = (i + DF) * W;
  1010.       for (j = 1; j < width; j += DF)
  1011.     for (k = j, l = DF - 1; k < j + DF - 1; k++, l--)
  1012.       {
  1013.         a[0][3] = a[1][3] = a[2][3] = a[3][3] = 1.0;
  1014.         a[0][0] = (gdouble) (k - 1);
  1015.         a[1][0] = (gdouble) (k + l);
  1016.         a[2][0] = (gdouble) (j - 1);
  1017.         a[3][0] = (gdouble) (j + DF - 1);
  1018.         a[0][1] = a[1][1] = (gdouble) i;
  1019.         a[2][1] = a[3][1] = (gdouble) (i + DF);
  1020.         a[0][2] = a[0][0] * a[0][1];
  1021.         a[1][2] = a[1][0] * a[1][1];
  1022.         a[2][2] = a[2][0] * a[2][1];
  1023.         a[3][2] = a[3][0] * a[3][1];
  1024.         b[0] = temp[(k - 1) + addr];
  1025.         b[1] = temp[(k + l) + addr];
  1026.         b[2] = temp[(j - 1) + addr2];
  1027.         b[3] = temp[(j + DF - 1) + addr2];
  1028.         dgefac (a, b);
  1029.         temp[k + addr] = (b[0] * (gdouble) k + b[1] * (gdouble) i +
  1030.                   b[2] * (gdouble) (i * k) + b[3]);
  1031.       }
  1032.     }
  1033.   /*
  1034.    * Column bilinear interpolation 
  1035.    */
  1036.   for (i = 1; i < height; i += DF)
  1037.     for (j = 0; j < width; j++)
  1038.       for (k = i, l = DF - 1; k < i + DF - 1; k++, l--)
  1039.     {
  1040.       a[0][3] = a[1][3] = a[2][3] = a[3][3] = 1.0;
  1041.       a[0][0] = a[2][0] = (gdouble) j;
  1042.       a[1][0] = a[3][0] = (gdouble) (j + 1);
  1043.       a[0][1] = a[1][1] = (gdouble) (k - 1);
  1044.       a[2][1] = a[3][1] = (gdouble) (k + l);
  1045.       a[0][2] = a[0][0] * a[0][1];
  1046.       a[1][2] = a[1][0] * a[1][1];
  1047.       a[2][2] = a[2][0] * a[2][1];
  1048.       a[3][2] = a[3][0] * a[3][1];
  1049.       b[0] = temp[j + (k - 1) * W];
  1050.       b[1] = temp[(j + 1) + (k - 1) * W];
  1051.       b[2] = temp[j + (k + l) * W];
  1052.       b[3] = temp[(j + 1) + (k + l) * W];
  1053.       dgefac (a, b);
  1054.       temp[j + k * W] = (b[0] * (gdouble) j + b[1] * (gdouble) k +
  1055.                  b[2] * (gdouble) (k * j) + b[3]);
  1056.     }
  1057.   for (i = 0; i < height; i++)
  1058.     {
  1059.       addr = i * W;
  1060.       for (j = 0; j < width; j++)
  1061.     dest[j + i * width] = temp[j + addr];
  1062.     }
  1063.   g_free (temp);
  1064. }
  1065.