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 FUNCTION (
00014 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 ) { register short j, k, nlinrem, npixrem;
00029     char   msg[SLEN];
00030     double x, y, *elin=0, *epix=0;
00031 
00032     if (TIMES) TIMING(T_FUNCTION);
00033     if (NAMES) {
00034         MESSAGE('I',"");
00035         MESSAGE('I',"FUNCTION");
00036         MESSAGE('I',"");
00037         sprintf(msg," Image size: lines:           %d",nlin);
00038         MESSAGE('I',msg);
00039         sprintf(msg,"             pixels:          %d",npix);
00040         MESSAGE('I',msg);
00041         sprintf(msg," Function radius: lines:    %12.4e",yradius);
00042         MESSAGE('I',msg);
00043         sprintf(msg,"                  pixels:   %12.4e",xradius);
00044         MESSAGE('I',msg);
00045         if (option == CONE) {
00046             MESSAGE('I'," Cone function");
00047         } else if (option == NEGEXP) {
00048             MESSAGE('I'," Negative exponential function");
00049         } else if (option == GAUSS) {
00050             MESSAGE('I'," Gaussian function");
00051         } else if (option == BOX) {
00052             MESSAGE('I'," Box function");
00053         } else if (option == CYLINDER) {
00054             MESSAGE('I'," Cylinder function");
00055         } else if (option == HANNING) {
00056             MESSAGE('I'," Hanning function");
00057         } else if (option == HAMMING) {
00058             MESSAGE('I'," Hanning function");
00059         } else if (option == KAISER) {
00060             MESSAGE('I'," Kaiser function");
00061         }
00062         MESSAGE('I'," ...............");
00063     }
00064 
00065     /* check input */
00066     if (nlin <= 0) {
00067         MESSAGE('E'," Number of lines must be greater than zero.");
00068         goto the_end;
00069     } else if (npix <= 0) {
00070         MESSAGE('E'," Number of pixels/line must be greater than zero.");
00071         goto the_end;
00072     } else if (yradius <= 0.0) {
00073         MESSAGE('E'," Radius in lines of test function must be greater than zero.");
00074         goto the_end;
00075     } else if (xradius <= 0.0) {
00076         MESSAGE('E'," Radius in pixels/line of test function must be greater than zero.");
00077         goto the_end;
00078     } else if (option != CONE  &&  option != NEGEXP  &&  option != GAUSS  &&  option != BOX  &&  option != CYLINDER && option != HANNING && option != HAMMING && option != KAISER) {
00079         MESSAGE('E'," Option must be one of CONE, NEGEXP, GAUSS, BOX, CYLINDER, HANNING, HAMMING or KAISER.");
00080         goto the_end;
00081     }
00082 
00083     /* allocate memory for temporary storage */
00084     if (!(elin = (double *)malloc(((nlin/2)+(nlin%2))*sizeof(double)))) {
00085         MESSAGE('E'," Memory request failed.");
00086         goto the_end;
00087     }
00088     if (!(epix = (double *)malloc(((npix/2)+(npix%2))*sizeof(double)))) {
00089         MESSAGE('E'," Memory request failed.");
00090         goto the_end;
00091     }
00092 
00093     /* create image of appropriate size */
00094     if (!CHECKIMG(*out)) GETMEM(1,nlin,npix,out);
00095     if (!*out) goto the_end;
00096 
00097     /* even/odd size parameter */
00098     nlinrem = nlin%2;
00099     npixrem = npix%2;
00100 
00101     /* set up the function image based on the option */
00102     switch (option) {
00103     case CONE:
00104         /* symmetric quadrants */
00105         for (j=1-nlinrem; j<nlin/2; j++) {
00106             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00107             for (k=1-npixrem; k<npix/2; k++) {
00108                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00109                 (*out)->data[0][j][k] =
00110                 (*out)->data[0][j][npix-k-npixrem] =
00111                 (*out)->data[0][nlin-j-nlinrem][k] =
00112                 (*out)->data[0][nlin-j-nlinrem][npix-k-npixrem] =
00113                 (PIXEL)max(0.0,(double)(SMAX-SMIN)*(1.0-sqrt(x+y)+(double)SMIN));
00114             }
00115         }
00116         /* zero coordinates */
00117         for (j=0; j<nlin; j++) {
00118             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00119             (*out)->data[0][j][npix/2] = (PIXEL)max(0.0,(double)(SMAX-SMIN)*(1.0-sqrt(y)+(double)SMIN));
00120         }
00121         for (k=0; k<npix; k++) {
00122             x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00123             (*out)->data[0][nlin/2][k] = (PIXEL)max(0.0,(double)(SMAX-SMIN)*(1.0-sqrt(x)+(double)SMIN));
00124         }
00125         /* first column (if npix is even) */
00126         if (npixrem == 0) {
00127             x = ((double)(npix/2)*(double)(npix/2))/(xradius*xradius);
00128             for (j=0; j<nlin; j++) {
00129                 y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00130                 (*out)->data[0][j][0] = (PIXEL)max(0.0,(double)(SMAX-SMIN)*(1.0-sqrt(x+y)+(double)SMIN));
00131             }
00132         }
00133         /* first row (if nlin is even) */
00134         if (nlinrem == 0) {
00135             y = ((double)(nlin/2)*(double)(nlin/2))/(yradius*yradius);
00136             for (k=0; k<npix; k++) {
00137                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00138                 (*out)->data[0][0][k] = (PIXEL)max(0.0,(double)(SMAX-SMIN)*(1.0-sqrt(x+y)+(double)SMIN));
00139             }
00140         }
00141         break;
00142     case NEGEXP:
00143         /* symmetric quadrants */
00144         for (j=1-nlinrem; j<nlin/2; j++) {
00145             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00146             for (k=1-npixrem; k<npix/2; k++) {
00147                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00148                 (*out)->data[0][j][k] =
00149                     (*out)->data[0][j][npix-k-npixrem] =
00150                 (*out)->data[0][nlin-j-nlinrem][k] =
00151                 (*out)->data[0][nlin-j-nlinrem][npix-k-npixrem] =
00152                 (PIXEL)((double)(SMAX-SMIN)*exp(-sqrt(x+y))+(double)SMIN);
00153             }
00154         }
00155         /* zero coordinates */
00156         for (j=0; j<nlin; j++) {
00157             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00158             (*out)->data[0][j][npix/2] = (PIXEL)((double)(SMAX-SMIN)*exp(-sqrt(y))+(double)SMIN);
00159         }
00160         for (k=0; k<npix; k++) {
00161             x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00162             (*out)->data[0][nlin/2][k] = (PIXEL)((double)(SMAX-SMIN)*exp(-sqrt(x))+(double)SMIN);
00163         }
00164         /* first column (if npix is even) */
00165         if (npixrem == 0) {
00166             x = ((double)(npix/2)*(double)(npix/2))/(xradius*xradius);
00167             for (j=0; j<nlin; j++) {
00168                 y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00169                 (*out)->data[0][j][0] = (PIXEL)((double)(SMAX-SMIN)*exp(-sqrt(x+y))+(double)SMIN);
00170             }
00171         }
00172         /* first row (if nlin is even) */
00173         if (nlinrem == 0) {
00174             y = ((double)(nlin/2)*(double)(nlin/2))/(yradius*yradius);
00175             for (k=0; k<npix; k++) {
00176                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00177                 (*out)->data[0][0][k] = (PIXEL)((double)(SMAX-SMIN)*exp(-sqrt(x+y))+(double)SMIN);
00178             }
00179         }
00180         break;
00181     case GAUSS:
00182         /* store one dimensional components */
00183         for (j=1-nlinrem; j<nlin/2; j++) {
00184             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00185             elin[j] = exp(-y);
00186         }
00187         for (k=1-npixrem; k<npix/2; k++) {
00188             x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00189             epix[k] = exp(-x);
00190         }
00191         /* symmetric quadrants */
00192         for (j=1-nlinrem; j<nlin/2; j++) {
00193             for (k=1-npixrem; k<npix/2; k++) {
00194                 (*out)->data[0][j][k] =
00195                 (*out)->data[0][j][npix-k-npixrem] =
00196                 (*out)->data[0][nlin-j-nlinrem][k] =
00197                 (*out)->data[0][nlin-j-nlinrem][npix-k-npixrem] =
00198                 (PIXEL)((double)(SMAX-SMIN)*(epix[k]*elin[j])+(double)SMIN);
00199             }
00200         }
00201         /* zero coordinates */
00202         for (j=0; j<nlin; j++) {
00203             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00204             (*out)->data[0][j][npix/2] = (PIXEL)((double)(SMAX-SMIN)*exp(-y)+(double)SMIN);
00205         }
00206         for (k=0; k<npix; k++) {
00207             x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00208             (*out)->data[0][nlin/2][k] = (PIXEL)((double)(SMAX-SMIN)*exp(-x)+(double)SMIN);
00209         }
00210         /* first column (if npix is even) */
00211         if (npixrem == 0) {
00212             x = ((double)(npix/2)*(double)(npix/2))/(xradius*xradius);
00213             for (j=0; j<nlin; j++) {
00214                 y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00215                 (*out)->data[0][j][0] = (PIXEL)((double)(SMAX-SMIN)*exp(-x-y)+(double)SMIN);
00216             }
00217         }
00218         /* first row (if nlin is even) */
00219         if (nlinrem == 0) {
00220             y = ((double)(nlin/2)*(double)(nlin/2))/(yradius*yradius);
00221             for (k=0; k<npix; k++) {
00222                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00223                 (*out)->data[0][0][k] = (PIXEL)((double)(SMAX-SMIN)*exp(-x-y)+(double)SMIN);
00224             }
00225         }
00226         break;
00227     case BOX:
00228         /* symmetric quadrants */
00229         for (j=1-nlinrem; j<nlin/2; j++) {
00230             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00231             for (k=1-npixrem; k<npix/2; k++) {
00232                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00233                 (*out)->data[0][j][k] =
00234                 (*out)->data[0][j][npix-k-npixrem] =
00235                 (*out)->data[0][nlin-j-nlinrem][k] =
00236                 (*out)->data[0][nlin-j-nlinrem][npix-k-npixrem] =
00237                 (PIXEL)(x <= 1.0  &&  y <= 1.0  ?  SMAX : SMIN);
00238             }
00239         }
00240         /* zero coordinates */
00241         for (j=0; j<nlin; j++) {
00242             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00243             (*out)->data[0][j][npix/2] = (PIXEL)(y <= 1.0  ?  SMAX : SMIN);
00244         }
00245         for (k=0; k<npix; k++) {
00246             x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00247             (*out)->data[0][nlin/2][k] = (PIXEL)(x <= 1.0  ?  SMAX : SMIN);
00248         }
00249         /* first column (if npix is even) */
00250         if (npixrem == 0) {
00251             x = ((double)(npix/2)*(double)(npix/2))/(xradius*xradius);
00252             for (j=0; j<nlin; j++) {
00253                 y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00254                 (*out)->data[0][j][0] = (PIXEL)(x <= 1.0  &&  y <= 1.0  ?  SMAX : SMIN);
00255             }
00256         }
00257         /* first row (if nlin is even) */
00258         if (nlinrem == 0) {
00259             y = ((double)(nlin/2)*(double)(nlin/2))/(yradius*yradius);
00260             for (k=0; k<npix; k++) {
00261                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00262                 (*out)->data[0][0][k] = (PIXEL)(x <= 1.0  &&  y <= 1.0  ?  SMAX : SMIN);
00263             }
00264         }
00265         break;
00266     case CYLINDER:
00267         /* symmetric quadrants */
00268         for (j=1-nlinrem; j<nlin/2; j++) {
00269             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00270             for (k=1-npixrem; k<npix/2; k++) {
00271                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00272                 (*out)->data[0][j][k] =
00273                 (*out)->data[0][j][npix-k-npixrem] =
00274                 (*out)->data[0][nlin-j-nlinrem][k] =
00275                 (*out)->data[0][nlin-j-nlinrem][npix-k-npixrem] =
00276                 (PIXEL)(sqrt(x+y) <= 1.0  ?  SMAX : SMIN);
00277             }
00278         }
00279         /* zero coordinates */
00280         for (j=0; j<nlin; j++) {
00281             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00282             (*out)->data[0][j][npix/2] = (PIXEL)(sqrt(y) <= 1.0  ?  SMAX : SMIN);
00283         }
00284         for (k=0; k<npix; k++) {
00285             x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00286             (*out)->data[0][nlin/2][k] = (PIXEL)(sqrt(x) <= 1.0  ?  SMAX : SMIN);
00287         }
00288         /* first column (if npix is even) */
00289         if (npixrem == 0) {
00290             x = ((double)(npix/2)*(double)(npix/2))/(xradius*xradius);
00291             for (j=0; j<nlin; j++) {
00292                 y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00293                 (*out)->data[0][j][0] = (PIXEL)(sqrt(x+y) <= 1.0  ?  SMAX : SMIN);
00294             }
00295         }
00296         /* first row (if nlin is even) */
00297         if (nlinrem == 0) {
00298             y = ((double)(nlin/2)*(double)(nlin/2))/(yradius*yradius);
00299             for (k=0; k<npix; k++) {
00300                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00301                 (*out)->data[0][0][k] = (PIXEL)(sqrt(x+y) <= 1.0  ?  SMAX : SMIN);
00302             }
00303         }
00304         break;
00305     case HANNING:
00306         /* symmetric quadrants */
00307         for (j=1-nlinrem; j<nlin/2; j++) {
00308             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00309             for (k=1-npixrem; k<npix/2; k++) {
00310                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00311                 (*out)->data[0][j][k] =
00312                 (*out)->data[0][j][npix-k-npixrem] =
00313                 (*out)->data[0][nlin-j-nlinrem][k] =
00314                 (*out)->data[0][nlin-j-nlinrem][npix-k-npixrem] =
00315                 (PIXEL)(sqrt(x+y) <= 1.0  ?  (SMAX-SMIN)*(0.50 + 0.50*cos(M_PI*sqrt(x+y))) + SMIN : SMIN);
00316             }
00317         }
00318         /* zero coordinates */
00319         for (j=0; j<nlin; j++) {
00320             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00321             (*out)->data[0][j][npix/2] = (PIXEL)(sqrt(y) <= 1.0  ?
00322                                                  (SMAX-SMIN)*(0.50 + 0.50*cos(M_PI*sqrt(y))) + SMIN : SMIN);
00323         }
00324         for (k=0; k<npix; k++) {
00325             x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00326             (*out)->data[0][nlin/2][k] = (PIXEL)(sqrt(x) <= 1.0  ?
00327                                                  (SMAX-SMIN)*(0.50 + 0.50*cos(M_PI*sqrt(x))) + SMIN : SMIN);
00328         }
00329         /* first column (if npix is even) */
00330         if (npixrem == 0) {
00331             x = ((double)(npix/2)*(double)(npix/2))/(xradius*xradius);
00332             for (j=0; j<nlin; j++) {
00333                 y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00334                 (*out)->data[0][j][0] = (PIXEL)(sqrt(x+y) <= 1.0  ?
00335                                                 (SMAX-SMIN)*(0.50 + 0.50*cos(M_PI*sqrt(x+y))) + SMIN : SMIN);
00336             }
00337         }
00338         /* first row (if nlin is even) */
00339         if (nlinrem == 0) {
00340             y = ((double)(nlin/2)*(double)(nlin/2))/(yradius*yradius);
00341             for (k=0; k<npix; k++) {
00342                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00343                 (*out)->data[0][0][k] = (PIXEL)(sqrt(x+y) <= 1.0  ?
00344                                                 (SMAX-SMIN)*(0.50 + 0.50*cos(M_PI*sqrt(x+y))) + SMIN : SMIN);
00345             }
00346         }
00347         break;
00348     case HAMMING:
00349         /* symmetric quadrants */
00350         for (j=1-nlinrem; j<nlin/2; j++) {
00351             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00352             for (k=1-npixrem; k<npix/2; k++) {
00353                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00354                 (*out)->data[0][j][k] =
00355                 (*out)->data[0][j][npix-k-npixrem] =
00356                 (*out)->data[0][nlin-j-nlinrem][k] =
00357                 (*out)->data[0][nlin-j-nlinrem][npix-k-npixrem] =
00358                 (PIXEL)(sqrt(x+y) <= 1.0  ?  (SMAX-SMIN)*(0.54 + 0.46*cos(M_PI*sqrt(x+y))) + SMIN : SMIN);
00359             }
00360         }
00361         /* zero coordinates */
00362         for (j=0; j<nlin; j++) {
00363             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00364             (*out)->data[0][j][npix/2] = (PIXEL)(sqrt(y) <= 1.0  ?
00365                                                  (SMAX-SMIN)*(0.54 + 0.46*cos(M_PI*sqrt(y))) + SMIN : SMIN);
00366         }
00367         for (k=0; k<npix; k++) {
00368             x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00369             (*out)->data[0][nlin/2][k] = (PIXEL)(sqrt(x) <= 1.0  ?
00370                                                  (SMAX-SMIN)*(0.54 + 0.46*cos(M_PI*sqrt(x))) + SMIN : SMIN);
00371         }
00372         /* first column (if npix is even) */
00373         if (npixrem == 0) {
00374             x = ((double)(npix/2)*(double)(npix/2))/(xradius*xradius);
00375             for (j=0; j<nlin; j++) {
00376                 y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00377                 (*out)->data[0][j][0] = (PIXEL)(sqrt(x+y) <= 1.0  ?
00378                                                 (SMAX-SMIN)*(0.54 + 0.46*cos(M_PI*sqrt(x+y))) + SMIN : SMIN);
00379             }
00380         }
00381         /* first row (if nlin is even) */
00382         if (nlinrem == 0) {
00383             y = ((double)(nlin/2)*(double)(nlin/2))/(yradius*yradius);
00384             for (k=0; k<npix; k++) {
00385                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00386                 (*out)->data[0][0][k] = (PIXEL)(sqrt(x+y) <= 1.0  ?
00387                                                 (SMAX-SMIN)*(0.54 + 0.46*cos(M_PI*sqrt(x+y))) + SMIN : SMIN);
00388             }
00389         }
00390         break;
00391     case KAISER:
00392         /* symmetric quadrants */
00393         for (j=1-nlinrem; j<nlin/2; j++) {
00394             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00395             for (k=1-npixrem; k<npix/2; k++) {
00396                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00397                 (*out)->data[0][j][k] =
00398                 (*out)->data[0][j][npix-k-npixrem] =
00399                 (*out)->data[0][nlin-j-nlinrem][k] =
00400                 (*out)->data[0][nlin-j-nlinrem][npix-k-npixrem] =
00401                 (PIXEL)(sqrt(x+y) <= 1.0  ?  (SMAX-SMIN)*(mBesselF(alpha*sqrt(1.0-(x+y)))/mBesselF(alpha)) + SMIN : SMIN);
00402             }
00403         }
00404         /* zero coordinates */
00405         for (j=0; j<nlin; j++) {
00406             y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00407             (*out)->data[0][j][npix/2] = (PIXEL)(sqrt(y) <= 1.0  ?
00408                                                  (SMAX-SMIN)*(mBesselF(alpha*sqrt(1.0-y))/mBesselF(alpha)) + SMIN : SMIN);
00409         }
00410         for (k=0; k<npix; k++) {
00411             x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00412             (*out)->data[0][nlin/2][k] = (PIXEL)(sqrt(x) <= 1.0  ?
00413                                                  (SMAX-SMIN)*(mBesselF(alpha*sqrt(1.0-x))/mBesselF(alpha)) + SMIN : SMIN);
00414         }
00415         /* first column (if npix is even) */
00416         if (npixrem == 0) {
00417             x = ((double)(npix/2)*(double)(npix/2))/(xradius*xradius);
00418             for (j=0; j<nlin; j++) {
00419                 y = ((double)(j-nlin/2)*(double)(j-nlin/2))/(yradius*yradius);
00420                 (*out)->data[0][j][0] = (PIXEL)(sqrt(x+y) <= 1.0  ?
00421                                                 (SMAX-SMIN)*(mBesselF(alpha*sqrt(1.0-(x+y)))/mBesselF(alpha)) + SMIN : SMIN);
00422             }
00423         }
00424         /* first row (if nlin is even) */
00425         if (nlinrem == 0) {
00426             y = ((double)(nlin/2)*(double)(nlin/2))/(yradius*yradius);
00427             for (k=0; k<npix; k++) {
00428                 x = ((double)(k-npix/2)*(double)(k-npix/2))/(xradius*xradius);
00429                 (*out)->data[0][0][k] = (PIXEL)(sqrt(x+y) <= 1.0  ?
00430                                                 (SMAX-SMIN)*(mBesselF(alpha*sqrt(1.0-(x+y)))/mBesselF(alpha)) + SMIN : SMIN);
00431             }
00432         }
00433         break;
00434     }
00435     (*out)->gmin = (*out)->data[0][0][0];
00436     (*out)->gmax = (*out)->data[0][nlin/2][npix/2];
00437 
00438     the_end:
00439     if (elin)  free(elin);
00440     if (epix)  free(epix);
00441     if (TIMES) TIMING(T_EXIT);
00442 }