function.c Source File

00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   This function creates an image as a two-dimensional test function.       */
00009 /*   The function maximum graylevel is set to SMAX.                           */
00010 /*                                                                            */
00011 /*----------------------------------------------------------------------------*/
00012 /*-Interface Information------------------------------------------------------*/
00013 void
00014 FUNCTION (short nlin,           /*  I   Number of lines in the output image.                  */
00015           short npix,           /*  I   Number of pixels/line in the output image.            */
00016           double yradius,       /*  I   Radius of the test function in y, in lines.           */
00017           double xradius,       /*  I   Radius of the test function in x, in pixels.          */
00018           double alpha,         /*  I   parameter for Kaiser window.                          */
00019           short option,         /*  I   Switch, indicating the type of the test function      */
00020           /*      ("r" is the distance of a pixel from the center):     */
00021           /*      CONE       -   SMAX * (1.0-r)                         */
00022           /*      NEGEXP     -   SMAX * exp(-r)                         */
00023           /*      GAUSS      -   SMAX * exp(-r**2)                      */
00024           /*      BOX        -   2*yradius lines by 2*xradius pixels    */
00025           /*      CYLINDER   -   r = 1.0                                */
00026           IMAGE ** out          /*  O   Address of a pointer to the output image.             */
00027 /*----------------------------------------------------------------------------*/
00028   )
00029 {
00030   register short j, k, nlinrem, npixrem;
00031   char msg[SLEN];
00032   double x, y, *elin = 0, *epix = 0;
00033 
00034   if (TIMES)
00035     TIMING (T_FUNCTION);
00036   if (NAMES)
00037     {
00038       MESSAGE ('I', "");
00039       MESSAGE ('I', "FUNCTION");
00040       MESSAGE ('I', "");
00041       sprintf (msg, " Image size: lines:           %d", nlin);
00042       MESSAGE ('I', msg);
00043       sprintf (msg, "             pixels:          %d", npix);
00044       MESSAGE ('I', msg);
00045       sprintf (msg, " Function radius: lines:    %12.4e", yradius);
00046       MESSAGE ('I', msg);
00047       sprintf (msg, "                  pixels:   %12.4e", xradius);
00048       MESSAGE ('I', msg);
00049       if (option == CONE)
00050         {
00051           MESSAGE ('I', " Cone function");
00052         }
00053       else if (option == NEGEXP)
00054         {
00055           MESSAGE ('I', " Negative exponential function");
00056         }
00057       else if (option == GAUSS)
00058         {
00059           MESSAGE ('I', " Gaussian function");
00060         }
00061       else if (option == BOX)
00062         {
00063           MESSAGE ('I', " Box function");
00064         }
00065       else if (option == CYLINDER)
00066         {
00067           MESSAGE ('I', " Cylinder function");
00068         }
00069       else if (option == HANNING)
00070         {
00071           MESSAGE ('I', " Hanning function");
00072         }
00073       else if (option == HAMMING)
00074         {
00075           MESSAGE ('I', " Hanning function");
00076         }
00077       else if (option == KAISER)
00078         {
00079           MESSAGE ('I', " Kaiser function");
00080         }
00081       MESSAGE ('I', " ...............");
00082     }
00083 
00084   /* check input */
00085   if (nlin <= 0)
00086     {
00087       MESSAGE ('E', " Number of lines must be greater than zero.");
00088       goto the_end;
00089     }
00090   else if (npix <= 0)
00091     {
00092       MESSAGE ('E', " Number of pixels/line must be greater than zero.");
00093       goto the_end;
00094     }
00095   else if (yradius <= 0.0)
00096     {
00097       MESSAGE ('E',
00098                " Radius in lines of test function must be greater than zero.");
00099       goto the_end;
00100     }
00101   else if (xradius <= 0.0)
00102     {
00103       MESSAGE ('E',
00104                " Radius in pixels/line of test function must be greater than zero.");
00105       goto the_end;
00106     }
00107   else if (option != CONE && option != NEGEXP && option != GAUSS
00108            && option != BOX && option != CYLINDER && option != HANNING
00109            && option != HAMMING && option != KAISER)
00110     {
00111       MESSAGE ('E',
00112                " Option must be one of CONE, NEGEXP, GAUSS, BOX, CYLINDER, HANNING, HAMMING or KAISER.");
00113       goto the_end;
00114     }
00115 
00116   /* allocate memory for temporary storage */
00117   if (!
00118       (elin =
00119        (double *) malloc (((nlin / 2) + (nlin % 2)) * sizeof (double))))
00120     {
00121       MESSAGE ('E', " Memory request failed.");
00122       goto the_end;
00123     }
00124   if (!
00125       (epix =
00126        (double *) malloc (((npix / 2) + (npix % 2)) * sizeof (double))))
00127     {
00128       MESSAGE ('E', " Memory request failed.");
00129       goto the_end;
00130     }
00131 
00132   /* create image of appropriate size */
00133   if (!CHECKIMG (*out))
00134     GETMEM (1, nlin, npix, out);
00135   if (!*out)
00136     goto the_end;
00137 
00138   /* even/odd size parameter */
00139   nlinrem = nlin % 2;
00140   npixrem = npix % 2;
00141 
00142   /* set up the function image based on the option */
00143   switch (option)
00144     {
00145     case CONE:
00146       /* symmetric quadrants */
00147       for (j = 1 - nlinrem; j < nlin / 2; j++)
00148         {
00149           y =
00150             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00151                                                                    yradius);
00152           for (k = 1 - npixrem; k < npix / 2; k++)
00153             {
00154               x =
00155                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00156                 (xradius * xradius);
00157               (*out)->data[0][j][k] = (*out)->data[0][j][npix - k - npixrem] =
00158                 (*out)->data[0][nlin - j - nlinrem][k] =
00159                 (*out)->data[0][nlin - j - nlinrem][npix - k - npixrem] =
00160                 (PIXEL) max (0.0,
00161                              (double) (SMAX - SMIN) * (1.0 - sqrt (x + y) +
00162                                                        (double) SMIN));
00163             }
00164         }
00165       /* zero coordinates */
00166       for (j = 0; j < nlin; j++)
00167         {
00168           y =
00169             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00170                                                                    yradius);
00171           (*out)->data[0][j][npix / 2] =
00172             (PIXEL) max (0.0,
00173                          (double) (SMAX - SMIN) * (1.0 - sqrt (y) +
00174                                                    (double) SMIN));
00175         }
00176       for (k = 0; k < npix; k++)
00177         {
00178           x =
00179             ((double) (k - npix / 2) * (double) (k - npix / 2)) / (xradius *
00180                                                                    xradius);
00181           (*out)->data[0][nlin / 2][k] =
00182             (PIXEL) max (0.0,
00183                          (double) (SMAX - SMIN) * (1.0 - sqrt (x) +
00184                                                    (double) SMIN));
00185         }
00186       /* first column (if npix is even) */
00187       if (npixrem == 0)
00188         {
00189           x =
00190             ((double) (npix / 2) * (double) (npix / 2)) / (xradius * xradius);
00191           for (j = 0; j < nlin; j++)
00192             {
00193               y =
00194                 ((double) (j - nlin / 2) * (double) (j - nlin / 2)) /
00195                 (yradius * yradius);
00196               (*out)->data[0][j][0] =
00197                 (PIXEL) max (0.0,
00198                              (double) (SMAX - SMIN) * (1.0 - sqrt (x + y) +
00199                                                        (double) SMIN));
00200             }
00201         }
00202       /* first row (if nlin is even) */
00203       if (nlinrem == 0)
00204         {
00205           y =
00206             ((double) (nlin / 2) * (double) (nlin / 2)) / (yradius * yradius);
00207           for (k = 0; k < npix; k++)
00208             {
00209               x =
00210                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00211                 (xradius * xradius);
00212               (*out)->data[0][0][k] =
00213                 (PIXEL) max (0.0,
00214                              (double) (SMAX - SMIN) * (1.0 - sqrt (x + y) +
00215                                                        (double) SMIN));
00216             }
00217         }
00218       break;
00219     case NEGEXP:
00220       /* symmetric quadrants */
00221       for (j = 1 - nlinrem; j < nlin / 2; j++)
00222         {
00223           y =
00224             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00225                                                                    yradius);
00226           for (k = 1 - npixrem; k < npix / 2; k++)
00227             {
00228               x =
00229                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00230                 (xradius * xradius);
00231               (*out)->data[0][j][k] = (*out)->data[0][j][npix - k - npixrem] =
00232                 (*out)->data[0][nlin - j - nlinrem][k] =
00233                 (*out)->data[0][nlin - j - nlinrem][npix - k - npixrem] =
00234                 (PIXEL) ((double) (SMAX - SMIN) * exp (-sqrt (x + y)) +
00235                          (double) SMIN);
00236             }
00237         }
00238       /* zero coordinates */
00239       for (j = 0; j < nlin; j++)
00240         {
00241           y =
00242             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00243                                                                    yradius);
00244           (*out)->data[0][j][npix / 2] =
00245             (PIXEL) ((double) (SMAX - SMIN) * exp (-sqrt (y)) +
00246                      (double) SMIN);
00247         }
00248       for (k = 0; k < npix; k++)
00249         {
00250           x =
00251             ((double) (k - npix / 2) * (double) (k - npix / 2)) / (xradius *
00252                                                                    xradius);
00253           (*out)->data[0][nlin / 2][k] =
00254             (PIXEL) ((double) (SMAX - SMIN) * exp (-sqrt (x)) +
00255                      (double) SMIN);
00256         }
00257       /* first column (if npix is even) */
00258       if (npixrem == 0)
00259         {
00260           x =
00261             ((double) (npix / 2) * (double) (npix / 2)) / (xradius * xradius);
00262           for (j = 0; j < nlin; j++)
00263             {
00264               y =
00265                 ((double) (j - nlin / 2) * (double) (j - nlin / 2)) /
00266                 (yradius * yradius);
00267               (*out)->data[0][j][0] =
00268                 (PIXEL) ((double) (SMAX - SMIN) * exp (-sqrt (x + y)) +
00269                          (double) SMIN);
00270             }
00271         }
00272       /* first row (if nlin is even) */
00273       if (nlinrem == 0)
00274         {
00275           y =
00276             ((double) (nlin / 2) * (double) (nlin / 2)) / (yradius * yradius);
00277           for (k = 0; k < npix; k++)
00278             {
00279               x =
00280                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00281                 (xradius * xradius);
00282               (*out)->data[0][0][k] =
00283                 (PIXEL) ((double) (SMAX - SMIN) * exp (-sqrt (x + y)) +
00284                          (double) SMIN);
00285             }
00286         }
00287       break;
00288     case GAUSS:
00289       /* store one dimensional components */
00290       for (j = 1 - nlinrem; j < nlin / 2; j++)
00291         {
00292           y =
00293             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00294                                                                    yradius);
00295           elin[j] = exp (-y);
00296         }
00297       for (k = 1 - npixrem; k < npix / 2; k++)
00298         {
00299           x =
00300             ((double) (k - npix / 2) * (double) (k - npix / 2)) / (xradius *
00301                                                                    xradius);
00302           epix[k] = exp (-x);
00303         }
00304       /* symmetric quadrants */
00305       for (j = 1 - nlinrem; j < nlin / 2; j++)
00306         {
00307           for (k = 1 - npixrem; k < npix / 2; k++)
00308             {
00309               (*out)->data[0][j][k] =
00310                 (*out)->data[0][j][npix - k - npixrem] =
00311                 (*out)->data[0][nlin - j - nlinrem][k] =
00312                 (*out)->data[0][nlin - j - nlinrem][npix - k - npixrem] =
00313                 (PIXEL) ((double) (SMAX - SMIN) * (epix[k] * elin[j]) +
00314                          (double) SMIN);
00315             }
00316         }
00317       /* zero coordinates */
00318       for (j = 0; j < nlin; j++)
00319         {
00320           y =
00321             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00322                                                                    yradius);
00323           (*out)->data[0][j][npix / 2] =
00324             (PIXEL) ((double) (SMAX - SMIN) * exp (-y) + (double) SMIN);
00325         }
00326       for (k = 0; k < npix; k++)
00327         {
00328           x =
00329             ((double) (k - npix / 2) * (double) (k - npix / 2)) / (xradius *
00330                                                                    xradius);
00331           (*out)->data[0][nlin / 2][k] =
00332             (PIXEL) ((double) (SMAX - SMIN) * exp (-x) + (double) SMIN);
00333         }
00334       /* first column (if npix is even) */
00335       if (npixrem == 0)
00336         {
00337           x =
00338             ((double) (npix / 2) * (double) (npix / 2)) / (xradius * xradius);
00339           for (j = 0; j < nlin; j++)
00340             {
00341               y =
00342                 ((double) (j - nlin / 2) * (double) (j - nlin / 2)) /
00343                 (yradius * yradius);
00344               (*out)->data[0][j][0] =
00345                 (PIXEL) ((double) (SMAX - SMIN) * exp (-x - y) +
00346                          (double) SMIN);
00347             }
00348         }
00349       /* first row (if nlin is even) */
00350       if (nlinrem == 0)
00351         {
00352           y =
00353             ((double) (nlin / 2) * (double) (nlin / 2)) / (yradius * yradius);
00354           for (k = 0; k < npix; k++)
00355             {
00356               x =
00357                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00358                 (xradius * xradius);
00359               (*out)->data[0][0][k] =
00360                 (PIXEL) ((double) (SMAX - SMIN) * exp (-x - y) +
00361                          (double) SMIN);
00362             }
00363         }
00364       break;
00365     case BOX:
00366       /* symmetric quadrants */
00367       for (j = 1 - nlinrem; j < nlin / 2; j++)
00368         {
00369           y =
00370             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00371                                                                    yradius);
00372           for (k = 1 - npixrem; k < npix / 2; k++)
00373             {
00374               x =
00375                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00376                 (xradius * xradius);
00377               (*out)->data[0][j][k] = (*out)->data[0][j][npix - k - npixrem] =
00378                 (*out)->data[0][nlin - j - nlinrem][k] =
00379                 (*out)->data[0][nlin - j - nlinrem][npix - k - npixrem] =
00380                 (PIXEL) (x <= 1.0 && y <= 1.0 ? SMAX : SMIN);
00381             }
00382         }
00383       /* zero coordinates */
00384       for (j = 0; j < nlin; j++)
00385         {
00386           y =
00387             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00388                                                                    yradius);
00389           (*out)->data[0][j][npix / 2] = (PIXEL) (y <= 1.0 ? SMAX : SMIN);
00390         }
00391       for (k = 0; k < npix; k++)
00392         {
00393           x =
00394             ((double) (k - npix / 2) * (double) (k - npix / 2)) / (xradius *
00395                                                                    xradius);
00396           (*out)->data[0][nlin / 2][k] = (PIXEL) (x <= 1.0 ? SMAX : SMIN);
00397         }
00398       /* first column (if npix is even) */
00399       if (npixrem == 0)
00400         {
00401           x =
00402             ((double) (npix / 2) * (double) (npix / 2)) / (xradius * xradius);
00403           for (j = 0; j < nlin; j++)
00404             {
00405               y =
00406                 ((double) (j - nlin / 2) * (double) (j - nlin / 2)) /
00407                 (yradius * yradius);
00408               (*out)->data[0][j][0] = (PIXEL) (x <= 1.0
00409                                                && y <= 1.0 ? SMAX : SMIN);
00410             }
00411         }
00412       /* first row (if nlin is even) */
00413       if (nlinrem == 0)
00414         {
00415           y =
00416             ((double) (nlin / 2) * (double) (nlin / 2)) / (yradius * yradius);
00417           for (k = 0; k < npix; k++)
00418             {
00419               x =
00420                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00421                 (xradius * xradius);
00422               (*out)->data[0][0][k] = (PIXEL) (x <= 1.0
00423                                                && y <= 1.0 ? SMAX : SMIN);
00424             }
00425         }
00426       break;
00427     case CYLINDER:
00428       /* symmetric quadrants */
00429       for (j = 1 - nlinrem; j < nlin / 2; j++)
00430         {
00431           y =
00432             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00433                                                                    yradius);
00434           for (k = 1 - npixrem; k < npix / 2; k++)
00435             {
00436               x =
00437                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00438                 (xradius * xradius);
00439               (*out)->data[0][j][k] = (*out)->data[0][j][npix - k - npixrem] =
00440                 (*out)->data[0][nlin - j - nlinrem][k] =
00441                 (*out)->data[0][nlin - j - nlinrem][npix - k - npixrem] =
00442                 (PIXEL) (sqrt (x + y) <= 1.0 ? SMAX : SMIN);
00443             }
00444         }
00445       /* zero coordinates */
00446       for (j = 0; j < nlin; j++)
00447         {
00448           y =
00449             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00450                                                                    yradius);
00451           (*out)->data[0][j][npix / 2] =
00452             (PIXEL) (sqrt (y) <= 1.0 ? SMAX : SMIN);
00453         }
00454       for (k = 0; k < npix; k++)
00455         {
00456           x =
00457             ((double) (k - npix / 2) * (double) (k - npix / 2)) / (xradius *
00458                                                                    xradius);
00459           (*out)->data[0][nlin / 2][k] =
00460             (PIXEL) (sqrt (x) <= 1.0 ? SMAX : SMIN);
00461         }
00462       /* first column (if npix is even) */
00463       if (npixrem == 0)
00464         {
00465           x =
00466             ((double) (npix / 2) * (double) (npix / 2)) / (xradius * xradius);
00467           for (j = 0; j < nlin; j++)
00468             {
00469               y =
00470                 ((double) (j - nlin / 2) * (double) (j - nlin / 2)) /
00471                 (yradius * yradius);
00472               (*out)->data[0][j][0] =
00473                 (PIXEL) (sqrt (x + y) <= 1.0 ? SMAX : SMIN);
00474             }
00475         }
00476       /* first row (if nlin is even) */
00477       if (nlinrem == 0)
00478         {
00479           y =
00480             ((double) (nlin / 2) * (double) (nlin / 2)) / (yradius * yradius);
00481           for (k = 0; k < npix; k++)
00482             {
00483               x =
00484                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00485                 (xradius * xradius);
00486               (*out)->data[0][0][k] =
00487                 (PIXEL) (sqrt (x + y) <= 1.0 ? SMAX : SMIN);
00488             }
00489         }
00490       break;
00491     case HANNING:
00492       /* symmetric quadrants */
00493       for (j = 1 - nlinrem; j < nlin / 2; j++)
00494         {
00495           y =
00496             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00497                                                                    yradius);
00498           for (k = 1 - npixrem; k < npix / 2; k++)
00499             {
00500               x =
00501                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00502                 (xradius * xradius);
00503               (*out)->data[0][j][k] = (*out)->data[0][j][npix - k - npixrem] =
00504                 (*out)->data[0][nlin - j - nlinrem][k] =
00505                 (*out)->data[0][nlin - j - nlinrem][npix - k - npixrem] =
00506                 (PIXEL) (sqrt (x + y) <=
00507                          1.0 ? (SMAX - SMIN) * (0.50 +
00508                                                 0.50 * cos (M_PI *
00509                                                             sqrt (x + y))) +
00510                          SMIN : SMIN);
00511             }
00512         }
00513       /* zero coordinates */
00514       for (j = 0; j < nlin; j++)
00515         {
00516           y =
00517             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00518                                                                    yradius);
00519           (*out)->data[0][j][npix / 2] =
00520             (PIXEL) (sqrt (y) <=
00521                      1.0 ? (SMAX - SMIN) * (0.50 +
00522                                             0.50 * cos (M_PI * sqrt (y))) +
00523                      SMIN : SMIN);
00524         }
00525       for (k = 0; k < npix; k++)
00526         {
00527           x =
00528             ((double) (k - npix / 2) * (double) (k - npix / 2)) / (xradius *
00529                                                                    xradius);
00530           (*out)->data[0][nlin / 2][k] =
00531             (PIXEL) (sqrt (x) <=
00532                      1.0 ? (SMAX - SMIN) * (0.50 +
00533                                             0.50 * cos (M_PI * sqrt (x))) +
00534                      SMIN : SMIN);
00535         }
00536       /* first column (if npix is even) */
00537       if (npixrem == 0)
00538         {
00539           x =
00540             ((double) (npix / 2) * (double) (npix / 2)) / (xradius * xradius);
00541           for (j = 0; j < nlin; j++)
00542             {
00543               y =
00544                 ((double) (j - nlin / 2) * (double) (j - nlin / 2)) /
00545                 (yradius * yradius);
00546               (*out)->data[0][j][0] =
00547                 (PIXEL) (sqrt (x + y) <=
00548                          1.0 ? (SMAX - SMIN) * (0.50 +
00549                                                 0.50 * cos (M_PI *
00550                                                             sqrt (x + y))) +
00551                          SMIN : SMIN);
00552             }
00553         }
00554       /* first row (if nlin is even) */
00555       if (nlinrem == 0)
00556         {
00557           y =
00558             ((double) (nlin / 2) * (double) (nlin / 2)) / (yradius * yradius);
00559           for (k = 0; k < npix; k++)
00560             {
00561               x =
00562                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00563                 (xradius * xradius);
00564               (*out)->data[0][0][k] =
00565                 (PIXEL) (sqrt (x + y) <=
00566                          1.0 ? (SMAX - SMIN) * (0.50 +
00567                                                 0.50 * cos (M_PI *
00568                                                             sqrt (x + y))) +
00569                          SMIN : SMIN);
00570             }
00571         }
00572       break;
00573     case HAMMING:
00574       /* symmetric quadrants */
00575       for (j = 1 - nlinrem; j < nlin / 2; j++)
00576         {
00577           y =
00578             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00579                                                                    yradius);
00580           for (k = 1 - npixrem; k < npix / 2; k++)
00581             {
00582               x =
00583                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00584                 (xradius * xradius);
00585               (*out)->data[0][j][k] = (*out)->data[0][j][npix - k - npixrem] =
00586                 (*out)->data[0][nlin - j - nlinrem][k] =
00587                 (*out)->data[0][nlin - j - nlinrem][npix - k - npixrem] =
00588                 (PIXEL) (sqrt (x + y) <=
00589                          1.0 ? (SMAX - SMIN) * (0.54 +
00590                                                 0.46 * cos (M_PI *
00591                                                             sqrt (x + y))) +
00592                          SMIN : SMIN);
00593             }
00594         }
00595       /* zero coordinates */
00596       for (j = 0; j < nlin; j++)
00597         {
00598           y =
00599             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00600                                                                    yradius);
00601           (*out)->data[0][j][npix / 2] =
00602             (PIXEL) (sqrt (y) <=
00603                      1.0 ? (SMAX - SMIN) * (0.54 +
00604                                             0.46 * cos (M_PI * sqrt (y))) +
00605                      SMIN : SMIN);
00606         }
00607       for (k = 0; k < npix; k++)
00608         {
00609           x =
00610             ((double) (k - npix / 2) * (double) (k - npix / 2)) / (xradius *
00611                                                                    xradius);
00612           (*out)->data[0][nlin / 2][k] =
00613             (PIXEL) (sqrt (x) <=
00614                      1.0 ? (SMAX - SMIN) * (0.54 +
00615                                             0.46 * cos (M_PI * sqrt (x))) +
00616                      SMIN : SMIN);
00617         }
00618       /* first column (if npix is even) */
00619       if (npixrem == 0)
00620         {
00621           x =
00622             ((double) (npix / 2) * (double) (npix / 2)) / (xradius * xradius);
00623           for (j = 0; j < nlin; j++)
00624             {
00625               y =
00626                 ((double) (j - nlin / 2) * (double) (j - nlin / 2)) /
00627                 (yradius * yradius);
00628               (*out)->data[0][j][0] =
00629                 (PIXEL) (sqrt (x + y) <=
00630                          1.0 ? (SMAX - SMIN) * (0.54 +
00631                                                 0.46 * cos (M_PI *
00632                                                             sqrt (x + y))) +
00633                          SMIN : SMIN);
00634             }
00635         }
00636       /* first row (if nlin is even) */
00637       if (nlinrem == 0)
00638         {
00639           y =
00640             ((double) (nlin / 2) * (double) (nlin / 2)) / (yradius * yradius);
00641           for (k = 0; k < npix; k++)
00642             {
00643               x =
00644                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00645                 (xradius * xradius);
00646               (*out)->data[0][0][k] =
00647                 (PIXEL) (sqrt (x + y) <=
00648                          1.0 ? (SMAX - SMIN) * (0.54 +
00649                                                 0.46 * cos (M_PI *
00650                                                             sqrt (x + y))) +
00651                          SMIN : SMIN);
00652             }
00653         }
00654       break;
00655     case KAISER:
00656       /* symmetric quadrants */
00657       for (j = 1 - nlinrem; j < nlin / 2; j++)
00658         {
00659           y =
00660             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00661                                                                    yradius);
00662           for (k = 1 - npixrem; k < npix / 2; k++)
00663             {
00664               x =
00665                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00666                 (xradius * xradius);
00667               (*out)->data[0][j][k] = (*out)->data[0][j][npix - k - npixrem] =
00668                 (*out)->data[0][nlin - j - nlinrem][k] =
00669                 (*out)->data[0][nlin - j - nlinrem][npix - k - npixrem] =
00670                 (PIXEL) (sqrt (x + y) <=
00671                          1.0 ? (SMAX -
00672                                 SMIN) * (mBesselF (alpha *
00673                                                    sqrt (1.0 -
00674                                                          (x +
00675                                                           y))) /
00676                                          mBesselF (alpha)) + SMIN : SMIN);
00677             }
00678         }
00679       /* zero coordinates */
00680       for (j = 0; j < nlin; j++)
00681         {
00682           y =
00683             ((double) (j - nlin / 2) * (double) (j - nlin / 2)) / (yradius *
00684                                                                    yradius);
00685           (*out)->data[0][j][npix / 2] =
00686             (PIXEL) (sqrt (y) <=
00687                      1.0 ? (SMAX -
00688                             SMIN) * (mBesselF (alpha * sqrt (1.0 - y)) /
00689                                      mBesselF (alpha)) + SMIN : SMIN);
00690         }
00691       for (k = 0; k < npix; k++)
00692         {
00693           x =
00694             ((double) (k - npix / 2) * (double) (k - npix / 2)) / (xradius *
00695                                                                    xradius);
00696           (*out)->data[0][nlin / 2][k] =
00697             (PIXEL) (sqrt (x) <=
00698                      1.0 ? (SMAX -
00699                             SMIN) * (mBesselF (alpha * sqrt (1.0 - x)) /
00700                                      mBesselF (alpha)) + SMIN : SMIN);
00701         }
00702       /* first column (if npix is even) */
00703       if (npixrem == 0)
00704         {
00705           x =
00706             ((double) (npix / 2) * (double) (npix / 2)) / (xradius * xradius);
00707           for (j = 0; j < nlin; j++)
00708             {
00709               y =
00710                 ((double) (j - nlin / 2) * (double) (j - nlin / 2)) /
00711                 (yradius * yradius);
00712               (*out)->data[0][j][0] =
00713                 (PIXEL) (sqrt (x + y) <=
00714                          1.0 ? (SMAX -
00715                                 SMIN) * (mBesselF (alpha *
00716                                                    sqrt (1.0 -
00717                                                          (x +
00718                                                           y))) /
00719                                          mBesselF (alpha)) + SMIN : SMIN);
00720             }
00721         }
00722       /* first row (if nlin is even) */
00723       if (nlinrem == 0)
00724         {
00725           y =
00726             ((double) (nlin / 2) * (double) (nlin / 2)) / (yradius * yradius);
00727           for (k = 0; k < npix; k++)
00728             {
00729               x =
00730                 ((double) (k - npix / 2) * (double) (k - npix / 2)) /
00731                 (xradius * xradius);
00732               (*out)->data[0][0][k] =
00733                 (PIXEL) (sqrt (x + y) <=
00734                          1.0 ? (SMAX -
00735                                 SMIN) * (mBesselF (alpha *
00736                                                    sqrt (1.0 -
00737                                                          (x +
00738                                                           y))) /
00739                                          mBesselF (alpha)) + SMIN : SMIN);
00740             }
00741         }
00742       break;
00743     }
00744   (*out)->gmin = (*out)->data[0][0][0];
00745   (*out)->gmax = (*out)->data[0][nlin / 2][npix / 2];
00746 
00747 the_end:
00748   if (elin)
00749     free (elin);
00750   if (epix)
00751     free (epix);
00752   if (TIMES)
00753     TIMING (T_EXIT);
00754 }