home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / mitsch75.zip / scheme-7_5_17-src.zip / scheme-7.5.17 / src / microcode / image.c < prev    next >
C/C++ Source or Header  |  1999-01-02  |  37KB  |  1,221 lines
  1. /* -*-C-*-
  2.  
  3. $Id: image.c,v 9.34 1999/01/02 06:11:34 cph Exp $
  4.  
  5. Copyright (c) 1987-1999 Massachusetts Institute of Technology
  6.  
  7. This program is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU General Public License as published by
  9. the Free Software Foundation; either version 2 of the License, or (at
  10. your option) any later version.
  11.  
  12. This program is distributed in the hope that it will be useful, but
  13. WITHOUT ANY WARRANTY; without even the implied warranty of
  14. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15. General Public License for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with this program; if not, write to the Free Software
  19. Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20. */
  21.  
  22. #include "scheme.h"
  23. #include "prims.h"
  24. #include "array.h"
  25. #include <math.h>
  26.  
  27. void
  28. arg_image (arg_number, nrows, ncols, array)
  29.      int arg_number;
  30.      long * nrows;
  31.      long * ncols;
  32.      REAL ** array;
  33. {
  34.   fast SCHEME_OBJECT argument = (ARG_REF (arg_number));
  35.   fast SCHEME_OBJECT rest;
  36.   fast SCHEME_OBJECT first;
  37.   fast SCHEME_OBJECT second;
  38.   fast SCHEME_OBJECT third;
  39.   if (! (PAIR_P (argument))) goto loser;
  40.   first = (PAIR_CAR (argument));
  41.   if (! (UNSIGNED_FIXNUM_P (first))) goto loser;
  42.   rest = (PAIR_CDR (argument));
  43.   if (! (PAIR_P (rest))) goto loser;
  44.   second = (PAIR_CAR (rest));
  45.   if (! (UNSIGNED_FIXNUM_P (second))) goto loser;
  46.   rest = (PAIR_CDR (rest));
  47.   if (! (PAIR_P (rest))) goto loser;
  48.   third = (PAIR_CAR (rest));
  49.   if (! (ARRAY_P (third))) goto loser;
  50.   if ((PAIR_CDR (rest)) != EMPTY_LIST) goto loser;
  51.   (*nrows) = (UNSIGNED_FIXNUM_TO_LONG (first));
  52.   (*ncols) = (UNSIGNED_FIXNUM_TO_LONG (second));
  53.   (*array) = (ARRAY_CONTENTS (third));
  54.   return;
  55.  loser:
  56.   error_bad_range_arg (arg_number);
  57.   /* NOTREACHED */
  58. }
  59.  
  60. #define MAKE_IMAGE(nrows, ncols, array)                    \
  61.   (cons ((LONG_TO_UNSIGNED_FIXNUM (nrows)),                \
  62.      (cons ((LONG_TO_UNSIGNED_FIXNUM (ncols)),            \
  63.         (cons ((array), EMPTY_LIST))))))
  64.  
  65. static int
  66. read_byte (fp)
  67.      fast FILE * fp;
  68. {
  69.   int result = (getc (fp));
  70.   if (ferror (fp))
  71.     error_external_return ();
  72.   return (result);
  73. }
  74.  
  75. static int
  76. read_word (fp)
  77.      fast FILE * fp;
  78. {
  79.   int result = (getw (fp));
  80.   if (ferror (fp))
  81.     error_external_return ();
  82.   return (result);
  83. }
  84.  
  85. static void
  86. write_word (fp, datum)
  87.      fast FILE * fp;
  88.      int datum;
  89. {
  90.   if ((putw (datum, fp)) != 0)
  91.     error_external_return ();
  92.   return;
  93. }
  94.  
  95. static int
  96. read_2bint (fp)
  97.      fast FILE * fp;
  98. {
  99.   int msd = (getc (fp));
  100.   if (ferror (fp))
  101.     error_external_return ();
  102.   {
  103.     int lsd = (getc (fp));
  104.     if (ferror (fp))
  105.       error_external_return ();
  106.     {
  107.       int result = ((msd << 8) | lsd);
  108.       return (((result & (1 << 15)) == 0) ? result : ((-1 << 16) | result));
  109.     }
  110.   }
  111. }
  112.  
  113. static void
  114. write_2bint (fp, datum)
  115.      fast FILE * fp;
  116.      int datum;
  117. {
  118.   if (((putc (((datum >> 8) & 0xFF), fp)) == EOF) ||
  119.       ((putc ((datum & 0xFF), fp)) == EOF))
  120.     error_external_return ();
  121.   return;
  122. }
  123.  
  124.  
  125. DEFINE_PRIMITIVE ("IMAGE-READ-ASCII", Prim_read_image_ascii, 1, 1, 0)
  126. {
  127.   fast FILE * fp;
  128.   long nrows, ncols;
  129.   PRIMITIVE_HEADER (1);
  130.   CHECK_ARG (1, STRING_P);
  131.   if ( (fp = fopen((STRING_ARG (1)), "r")) == NULL) 
  132.     error_bad_range_arg (1);
  133.   fscanf (fp, "%d %d \n", (&nrows), (&ncols));
  134.   if ((ferror (fp)) || ((ncols > 512) || (nrows > 512)))
  135.     { printf("read-image-ascii-file: problem with rows,cols \n");
  136.       error_bad_range_arg (1); }
  137.   {
  138.     fast long length = (nrows * ncols);
  139.     SCHEME_OBJECT array = (allocate_array (length));
  140.     fast REAL * scan = (ARRAY_CONTENTS (array));
  141.     while ((length--) > 0)
  142.       { long number;
  143.     fscanf (fp, "%d", (&number));
  144.     if (ferror (fp))
  145.       error_external_return ();
  146.     (*scan++) = ((REAL) number);
  147.       }
  148.     if ((fclose (fp)) != 0) error_external_return ();
  149.     PRIMITIVE_RETURN (MAKE_IMAGE (nrows, ncols, array));
  150.   }
  151. }
  152.  
  153. DEFINE_PRIMITIVE ("IMAGE-READ-2BINT", Prim_read_image_2bint, 1, 1, 0)
  154. {
  155.   FILE *fp, *fopen();
  156.   PRIMITIVE_HEADER (1);
  157.   CHECK_ARG (1, STRING_P);
  158.   if ( ( fp = (fopen((STRING_ARG (1)), "r")) ) == NULL) 
  159.     error_bad_range_arg (1);
  160.   {
  161.     int nrows = (read_word (fp));
  162.     int ncols = (read_word (fp));
  163.     fast long length = (nrows * ncols);
  164.     SCHEME_OBJECT array = (allocate_array (length));
  165.     fast REAL * scan = (ARRAY_CONTENTS (array));
  166.     while ((length--) > 0)
  167.       (*scan++) = ((REAL) (read_2bint (fp)));
  168.     if ((fclose (fp)) != 0)
  169.       error_external_return ();
  170.     PRIMITIVE_RETURN (MAKE_IMAGE (nrows, ncols, array));
  171.   }
  172. }
  173.  
  174. DEFINE_PRIMITIVE ("IMAGE-READ-CTSCAN", Prim_read_image_ctscan, 1, 1, 0)
  175. {
  176.   fast FILE * fp;
  177.   PRIMITIVE_HEADER (1);
  178.   CHECK_ARG (1, STRING_P);
  179.   fp = (fopen((STRING_ARG (1)), "r"));
  180.   if (fp == ((FILE *) 0))
  181.     error_bad_range_arg (1);
  182.   Primitive_GC_If_Needed (BYTES_TO_WORDS (512 * (sizeof (int))));
  183.   {
  184.     int nrows = 512;
  185.     int ncols = 512;
  186.     SCHEME_OBJECT array = (allocate_array (nrows * ncols));
  187.     REAL * Array = (ARRAY_CONTENTS (array));
  188.     fast int * Widths = ((int *) Free);
  189.     fast int i;
  190.     /* Discard header */
  191.     for (i = 0; (i < 2048); i += 1)
  192.       (void) read_byte (fp);
  193.     for (i = 0; (i < 512); i += 1)
  194.       (Widths [i]) = (read_2bint (fp));
  195.     for (i = 0; (i < (nrows * ncols)); i += 1)
  196.       (Array [i]) = 0;
  197.     for (i = 0; (i < 512); i += 1)
  198.       {
  199.     fast int array_index = ((i * 512) + (256 - (Widths [i])));
  200.     fast int m;
  201.     for (m = 0; (m < (2 * (Widths [i]))); m += 1)
  202.       (Array [array_index + m]) = ((REAL) (read_2bint (fp)));
  203.       }
  204.     /* CTSCAN images are upside down */
  205. #if (REAL_IS_DEFINED_DOUBLE != 0)
  206.     ALIGN_FLOAT (Free);
  207.     Free += 1;
  208. #endif
  209.     Primitive_GC_If_Needed (512 * REAL_SIZE);
  210.     Image_Mirror_Upside_Down (Array, nrows, ncols, ((REAL *) Free));
  211.     if ((fclose (fp)) != 0)
  212.       error_external_return ();
  213.     PRIMITIVE_RETURN (MAKE_IMAGE (nrows, ncols, array));
  214.   }
  215. }
  216.  
  217.  
  218. DEFINE_PRIMITIVE ("IMAGE-READ-CBIN", Prim_read_image_cbin, 1, 1, 0)
  219. {
  220.   fast FILE * fp;
  221.   PRIMITIVE_HEADER (1);
  222.   CHECK_ARG (1, STRING_P);
  223.   fp = (fopen ((STRING_ARG (1)), "r"));
  224.   if (fp == ((FILE *) 0))
  225.     error_bad_range_arg (1);
  226.   {
  227.     int nrows = (read_word (fp));
  228.     int ncols = (read_word (fp));
  229.     long length = (nrows * ncols);
  230.     SCHEME_OBJECT array = (allocate_array (length));
  231.     fast REAL * scan = (ARRAY_CONTENTS (array));
  232.     while ((length--) > 0)
  233.       (*scan++) = (read_word (fp));
  234.     if ((fclose (fp)) != 0)
  235.       error_external_return ();
  236.     PRIMITIVE_RETURN (MAKE_IMAGE (nrows, ncols, array));
  237.   }
  238. }
  239.  
  240.  
  241. Image_Mirror_Upside_Down (Array,nrows,ncols,Temp_Row)
  242.      REAL * Array;
  243.      long nrows;
  244.      long ncols;
  245.      REAL * Temp_Row;
  246. { int i;
  247.   REAL *M_row, *N_row;
  248.   for (i=0;i<(nrows/2);i++) {
  249.     M_row = Array + (i * ncols);
  250.     N_row = Array + (((nrows-1)-i) * ncols);
  251.     C_Array_Copy(N_row,    Temp_Row, ncols);
  252.     C_Array_Copy(M_row,    N_row,    ncols);
  253.     C_Array_Copy(Temp_Row, M_row,    ncols);
  254.   }
  255. }
  256.  
  257. /* The following does not work, to be fixed. */
  258. DEFINE_PRIMITIVE ("IMAGE-DOUBLE-TO-FLOAT!", Prim_image_double_to_float, 1, 1, 0)
  259. {
  260.   long Length;
  261.   long i,j;
  262.   long allocated_cells;
  263.   long nrows, ncols;
  264.   SCHEME_OBJECT array;
  265.   double *Array, *From_Here;
  266.   fast double temp_value_cell;
  267.   float *To_Here;
  268.   PRIMITIVE_HEADER (1);
  269.   arg_image (1, (&nrows), (&ncols), (&array));
  270.   Array = ((double *) (ARRAY_CONTENTS (array)));
  271.   From_Here = Array;
  272.   To_Here = ((float *) (Array));
  273.   Length = nrows * ncols;
  274.  
  275.   for (i=0;i<Length;i++) {
  276.     temp_value_cell = *From_Here;
  277.     From_Here++;
  278.     *To_Here = ((float) temp_value_cell);
  279.     To_Here++;
  280.   }
  281.   /* and now SIDE-EFFECT the ARRAY_HEADER */
  282.   SET_VECTOR_LENGTH (array, ((Length * FLOAT_SIZE) + 1));
  283.   PRIMITIVE_RETURN (UNSPECIFIC);
  284. }
  285.  
  286.  
  287. DEFINE_PRIMITIVE ("SUBIMAGE-COPY!", Prim_subimage_copy, 12, 12, 0)
  288. {
  289.   long r1, c1, r2, c2, at1r, at1c, at2r, at2c, mr, mc;
  290.   REAL * x;
  291.   REAL * y;
  292.   void subimage_copy ();
  293.   PRIMITIVE_HEADER (12);
  294.   CHECK_ARG (3, ARRAY_P);    /* image array 1 = source array */
  295.   CHECK_ARG (6, ARRAY_P);    /* image array 2 = destination array */
  296.   r1 = (arg_nonnegative_integer (1));
  297.   c1 = (arg_nonnegative_integer (2));
  298.   if ((r1 * c1) != (ARRAY_LENGTH (ARG_REF (3))))
  299.     error_bad_range_arg (3);
  300.   x = (ARRAY_CONTENTS (ARG_REF (3)));
  301.   r2 = arg_nonnegative_integer (4);
  302.   c2 = arg_nonnegative_integer (5);
  303.   if ((r2 * c2) != (ARRAY_LENGTH (ARG_REF (6))))
  304.     error_bad_range_arg (6);
  305.   y = (ARRAY_CONTENTS (ARG_REF (6)));
  306.   at1r = (arg_nonnegative_integer (7));
  307.   at1c = (arg_nonnegative_integer (8));
  308.   at2r = (arg_nonnegative_integer (9));
  309.   at2c = (arg_nonnegative_integer (10));
  310.   mr = (arg_nonnegative_integer (11));
  311.   mc = (arg_nonnegative_integer (12));
  312.   if (((at1r + mr) > r1) || ((at1c + mc) > c1))
  313.     error_bad_range_arg (7);
  314.   if (((at2r + mr) > r2) || ((at2c + mc) > c2))
  315.     error_bad_range_arg (9);
  316.   subimage_copy (x, y, r1, c1, r2, c2, at1r, at1c, at2r, at2c, mr, mc);
  317.   PRIMITIVE_RETURN (UNSPECIFIC);
  318. }
  319.  
  320. void
  321. subimage_copy (x,y, r1,c1,r2,c2, at1r,at1c,at2r,at2c, mr,mc)
  322.      REAL *x,*y;
  323.      long r1,c1,r2,c2, at1r,at1c,at2r,at2c, mr,mc;
  324. { long i,j;
  325.   REAL *xrow,*yrow;
  326.  
  327.   xrow = x + at1r*c1   + at1c;
  328.   yrow = y + at2r*c2   + at2c;    /*  A(i,j)--->Array[i*ncols+j]  */
  329.  
  330.   for (i=0; i<mr; i++) {
  331.     for (j=0; j<mc; j++)    yrow[j] = xrow[j];
  332.     xrow = xrow + c1;
  333.     yrow = yrow + c2;
  334.   }
  335. }
  336.  
  337. /* image-operation-2
  338.    groups together procedures     that use 2 image-arrays
  339.    (usually side-effecting the 2nd image, but not necessarily) */
  340.  
  341. DEFINE_PRIMITIVE ("IMAGE-OPERATION-2!", Prim_image_operation_2, 5,5, 0)
  342. { long rows, cols, nn, opcode;
  343.   REAL *x,*y;
  344.   void image_laplacian();
  345.   PRIMITIVE_HEADER (5);
  346.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  347.   CHECK_ARG (2, FIXNUM_P);    /* rows */
  348.   CHECK_ARG (3, FIXNUM_P);    /* cols */
  349.   CHECK_ARG (4, ARRAY_P);    /* image array 1 */
  350.   CHECK_ARG (5, ARRAY_P);    /* image array 2 */
  351.  
  352.   x = ARRAY_CONTENTS(ARG_REF(4));
  353.   y = ARRAY_CONTENTS(ARG_REF(5));
  354.   rows = arg_nonnegative_integer(2);
  355.   cols = arg_nonnegative_integer(3);
  356.   nn = rows*cols;
  357.   if (nn != ARRAY_LENGTH(ARG_REF(4)))   error_bad_range_arg(4);
  358.   if (nn != ARRAY_LENGTH(ARG_REF(5)))   error_bad_range_arg(5);
  359.  
  360.   opcode = arg_nonnegative_integer(1);
  361.  
  362.   if (opcode==1)
  363.     image_laplacian(x,y,rows,cols); /* result in y */
  364.   else if (opcode==2)
  365.     error_bad_range_arg(1);    /* illegal opcode */
  366.   else
  367.     error_bad_range_arg(1);    /* illegal opcode */
  368.  
  369.   PRIMITIVE_RETURN (UNSPECIFIC);
  370. }
  371.  
  372. /* Laplacian form [4,-1,-1,-1,-1]/4
  373.    A(i,j) --> Array[i*ncols + j]
  374.    With no knowledge outside boundary,
  375.    assume laplace(edge-point)=0.0 (no wrap-around, no artificial bndry) */
  376. void
  377. image_laplacian (x,y, nrows,ncols)
  378.      REAL *x, *y;
  379.      long nrows, ncols;
  380. { long i,j, nrows1, ncols1;
  381.   nrows1=nrows-1; ncols1=ncols-1;
  382.   /* no need todo anything for 1-point image */
  383.   if ((nrows<2)||(ncols<2))
  384.     return;
  385.   /* */
  386.   i=0;j=0;           y[i*ncols+j] = 0.0; /* NE corner */
  387.   i=0;j=ncols1;      y[i*ncols+j] = 0.0; /* NW corner */
  388.   i=nrows1;j=0;      y[i*ncols+j] = 0.0; /* SE corner */
  389.   i=nrows1;j=ncols1; y[i*ncols+j] = 0.0; /* SW corner */
  390.   i=0; for (j=1;j<ncols1;j++)       y[i*ncols+j] = 0.0;    /* NORTH row */
  391.   i=nrows1; for (j=1;j<ncols1;j++)  y[i*ncols+j] = 0.0;    /* SOUTH row */
  392.   j=0; for (i=1;i<nrows1;i++)       y[i*ncols+j] = 0.0;    /* EAST column */
  393.   j=ncols1; for (i=1;i<nrows1;i++)  y[i*ncols+j] = 0.0;    /* WEST column */
  394.   /* */
  395.   for (i=1;i<nrows1;i++)
  396.     for (j=1;j<ncols1;j++)
  397.       y[i*ncols+j] =
  398.     x[i*ncols+j] -
  399.       ((x[i*ncols+(j-1)] +
  400.         x[i*ncols+(j+1)] +
  401.         x[(i-1)*ncols+j] +
  402.         x[(i+1)*ncols+j])
  403.        / 4);
  404.   return;
  405. }
  406.  
  407. DEFINE_PRIMITIVE ("IMAGE-DOUBLE-BY-INTERPOLATION", Prim_image_double_by_interpolation, 1, 1, 0)
  408. { long nrows, ncols, Length;
  409.   SCHEME_OBJECT Parray;
  410.   REAL *Array, *To_Here;
  411.   SCHEME_OBJECT Result, array;
  412.   PRIMITIVE_HEADER (1);
  413.   arg_image (1, (&nrows), (&ncols), (&Parray));
  414.   Array = (ARRAY_CONTENTS (Parray));
  415.   array = (allocate_array (4 * nrows * ncols));
  416.   Result = (MAKE_IMAGE (nrows, ncols, array));
  417.  
  418.   Array = ARRAY_CONTENTS(Parray);
  419.   C_image_double_by_interpolation
  420.     (Array, (ARRAY_CONTENTS(array)), nrows, ncols);
  421.   PRIMITIVE_RETURN(Result);
  422. }
  423.  
  424. /* double image by linear interpolation.
  425.    ---new_array must be 4 times as long ---
  426.    A(i,j) --> Array[i*ncols + j]
  427.    magnification in a south-east direction
  428.    (i.e. replication of pixels in South-East corner) */
  429. C_image_double_by_interpolation (array, new_array, nrows, ncols)
  430.      REAL *array, *new_array;
  431.      long nrows, ncols;
  432. { long i,j, nrows1, ncols1, nrows2, ncols2;
  433.   nrows1=nrows-1; ncols1=ncols-1;
  434.   nrows2=2*nrows; ncols2=2*ncols;
  435.   /* no need todo anything for 1-point image */
  436.   if ((nrows<2)||(ncols<2)) return(1);
  437.   i=nrows1; for (j=0;j<ncols1;j++)     /* SOUTH row */
  438.   { new_array[(2*i)*ncols2+(2*j)]      = array[i*ncols+j];
  439.     new_array[(2*i+1)*ncols2+(2*j)]    = array[i*ncols+j];
  440.     new_array[(2*i)*ncols2+(2*j)+1] =
  441.       ((array[i*ncols+j]+array[i*ncols+j+1]) / 2);
  442.     new_array[(2*i+1)*ncols2+(2*j)+1]  = new_array[(2*i)*ncols2+(2*j)+1];
  443.   }
  444.   j=ncols1; for (i=0;i<nrows1;i++)      /* WEST column */
  445.   { new_array[(2*i)*ncols2+(2*j)]      = array[i*ncols+j];
  446.     new_array[(2*i)*ncols2+(2*j)+1]    = array[i*ncols+j];
  447.     new_array[(2*i+1)*ncols2+(2*j)] =
  448.       ((array[i*ncols+j]+array[(i+1)*ncols+j]) / 2);
  449.     new_array[(2*i+1)*ncols2+(2*j)+1]  = new_array[(2*i+1)*ncols2+(2*j)];
  450.   }
  451.   i=nrows1;j=ncols1; {                  /* SW corner */
  452.     new_array[(2*i)*ncols2+(2*j)]     =  array[i*ncols+j];
  453.     new_array[(2*i)*ncols2+(2*j)+1]   =  array[i*ncols+j];
  454.     new_array[(2*i+1)*ncols2+(2*j)]   =  array[i*ncols+j];
  455.     new_array[(2*i+1)*ncols2+(2*j)+1] =  array[i*ncols+j];
  456.   }
  457.   /* */
  458.   for (i=0;i<nrows1;i++)
  459.     for (j=0;j<ncols1;j++) {                        /* interior of image */
  460.       new_array[(2*i)*ncols2+(2*j)] =  array[i*ncols+j];
  461.       new_array[(2*i)*ncols2+(2*j)+1] =
  462.     ((array[i*ncols+j]+array[i*ncols+j+1]) / 2);
  463.       new_array[(2*i+1)*ncols2+(2*j)] =
  464.     ((array[i*ncols+j]+array[(i+1)*ncols+j]) / 2);
  465.       new_array[(2*i+1)*ncols2+(2*j)+1] =
  466.     ((array[i*ncols+j] +
  467.       array[i*ncols+j+1] +
  468.       array[(i+1)*ncols+j] +
  469.       array[(i+1)*ncols+j+1])
  470.      / 4);
  471.     }
  472. }
  473.  
  474. DEFINE_PRIMITIVE ("IMAGE-MAKE-RING", Prim_image_make_ring, 4, 4, 0)
  475. {
  476.   PRIMITIVE_HEADER (4);
  477.   {
  478.     long nrows = (arg_nonnegative_integer (1));
  479.     long ncols = (arg_nonnegative_integer (2));
  480.     long Length = (nrows * ncols);
  481.     long high_cycle =
  482.       (arg_index_integer (4, ((min ((nrows / 2), (ncols / 2))) + 1)));
  483.     long low_cycle = (arg_index_integer (3, (high_cycle + 1)));
  484.     SCHEME_OBJECT Ring_Array_Result = (allocate_array (Length));
  485.     SCHEME_OBJECT Result = (MAKE_IMAGE (nrows, ncols, Ring_Array_Result));
  486.     REAL * Ring_Array = (ARRAY_CONTENTS (Ring_Array_Result));
  487.     C_Image_Make_Ring (Ring_Array, nrows, ncols, low_cycle, high_cycle);
  488.     PRIMITIVE_RETURN (Result);
  489.   }
  490. }
  491.  
  492. C_Image_Make_Ring (Ring_Array, nrows, ncols, low_cycle, high_cycle)
  493.      REAL *Ring_Array;
  494.      long nrows, ncols, low_cycle, high_cycle;
  495. { long Square_LC=low_cycle*low_cycle, Square_HC=high_cycle*high_cycle;
  496.   long i, j, m, n, radial_cycle;
  497.   long nrows2=nrows/2, ncols2=ncols/2;
  498.   for (i=0; i<nrows; i++) {
  499.     for (j=0; j<ncols; j++) {
  500.       m = ((i<nrows2) ? i : (nrows-i));
  501.       n = ((j<ncols2) ? j : (ncols-j));
  502.       radial_cycle = (m*m)+(n*n);
  503.       if ( (radial_cycle<Square_LC) || (radial_cycle>Square_HC))
  504.     Ring_Array[i*ncols+j] = 0;
  505.       else Ring_Array[i*ncols+j] = 1;
  506.     }}
  507. }
  508.  
  509. /* Periodic-shift without side-effects for code simplicity. */
  510.  
  511. DEFINE_PRIMITIVE ("IMAGE-PERIODIC-SHIFT", Prim_image_periodic_shift, 3, 3, 0)
  512. {
  513.   long nrows;
  514.   long ncols;
  515.   REAL * Array;
  516.   PRIMITIVE_HEADER (3);
  517.   {
  518.     SCHEME_OBJECT Parray;
  519.     arg_image (1, (&nrows), (&ncols), (&Parray));
  520.     Array = (ARRAY_CONTENTS (Parray));
  521.   }
  522.   {
  523.     long ver_shift = ((arg_integer (2)) % nrows);
  524.     long hor_shift = ((arg_integer (3)) % ncols);
  525.     SCHEME_OBJECT array = (allocate_array (nrows * ncols));
  526.     SCHEME_OBJECT Result = (MAKE_IMAGE (nrows, ncols, array));
  527.     C_Image_Periodic_Shift
  528.       (Array,
  529.        (ARRAY_CONTENTS (array)),
  530.        nrows,
  531.        ncols,
  532.        ver_shift,
  533.        hor_shift);
  534.     PRIMITIVE_RETURN (Result);
  535.   }
  536. }
  537.  
  538. /* ASSUMES ((hor_shift < nrows) && (ver_shift < ncols)) */
  539.  
  540. C_Image_Periodic_Shift(Array, New_Array, nrows, ncols, ver_shift, hor_shift)
  541.      REAL *Array, *New_Array; long nrows, ncols, hor_shift, ver_shift;
  542. { long i, j, ver_index, hor_index;
  543.   REAL *To_Here;
  544.   To_Here = New_Array;
  545.   for (i=0;i<nrows;i++) {
  546.     for (j=0;j<ncols;j++) {
  547.       ver_index = (i+ver_shift) % nrows;
  548.       if (ver_index<0) ver_index = nrows+ver_index; /* wrapping around */
  549.       hor_index = (j+hor_shift) % ncols;
  550.       if (hor_index<0) hor_index = ncols+hor_index;
  551.       *To_Here++ = Array[ver_index*ncols + hor_index];
  552.     }}
  553. }
  554.  
  555. /* Rotations and stuff */
  556.  
  557. DEFINE_PRIMITIVE ("IMAGE-TRANSPOSE!", Prim_image_transpose, 4,4, 0)
  558. { long rows, cols;
  559.   REAL *x, *y;
  560.   PRIMITIVE_HEADER (4);
  561.   CHECK_ARG (1, FIXNUM_P);    /* rows */
  562.   CHECK_ARG (2, FIXNUM_P);    /* cols */
  563.   CHECK_ARG (3, ARRAY_P);    /* image array 1 */
  564.   CHECK_ARG (4, ARRAY_P);    /* image array 2,  empty for rows=cols */
  565.   
  566.   rows = arg_nonnegative_integer(1);
  567.   cols = arg_nonnegative_integer(2);
  568.   x = (ARRAY_CONTENTS (ARG_REF(3)));
  569.   y = (ARRAY_CONTENTS (ARG_REF(4)));
  570.   
  571.   if (rows==cols)        /* square image ==> ignore argument 4  */
  572.     Image_Fast_Transpose (x, rows);
  573.   else
  574.     Image_Transpose (x, y, rows, cols);
  575.   
  576.   PRIMITIVE_RETURN (UNSPECIFIC);
  577. }
  578.  
  579. DEFINE_PRIMITIVE ("IMAGE-ROTATE-90CLW!", Prim_image_rotate_90clw, 1, 1, 0)
  580. {
  581.   long nrows;
  582.   long ncols;
  583.   REAL * Array;
  584.   long Length;
  585.   PRIMITIVE_HEADER (1);
  586.   {
  587.     SCHEME_OBJECT Parray;
  588.     arg_image (1, (&nrows), (&ncols), (&Parray));
  589.     Array = (ARRAY_CONTENTS (Parray));
  590.   }
  591.   Length = (nrows * ncols);
  592. #if (REAL_IS_DEFINED_DOUBLE != 0)
  593.     ALIGN_FLOAT (Free);
  594.     Free += 1;
  595. #endif
  596.   Primitive_GC_If_Needed (Length * REAL_SIZE);
  597.   {
  598.     REAL * Temp_Array = ((REAL *) Free);
  599.     Image_Rotate_90clw (Array, Temp_Array, nrows, ncols);
  600.     C_Array_Copy (Temp_Array, Array, Length);
  601.   }
  602.   {
  603.     SCHEME_OBJECT argument = (ARG_REF (1));
  604.     SET_PAIR_CAR (argument, (LONG_TO_UNSIGNED_FIXNUM (ncols)));
  605.     SET_PAIR_CAR ((PAIR_CDR (argument)), (LONG_TO_UNSIGNED_FIXNUM (nrows)));
  606.     PRIMITIVE_RETURN (argument);
  607.   }
  608. }
  609.  
  610. DEFINE_PRIMITIVE ("IMAGE-ROTATE-90CCLW!", Prim_image_rotate_90cclw, 1, 1, 0)
  611. {
  612.   long nrows;
  613.   long ncols;
  614.   REAL * Array;
  615.   long Length;
  616.   PRIMITIVE_HEADER (1);
  617.   {
  618.     SCHEME_OBJECT Parray;
  619.     arg_image (1, (&nrows), (&ncols), (&Parray));
  620.     Array = (ARRAY_CONTENTS (Parray));
  621.   }
  622.   Length = (nrows * ncols);
  623. #if (REAL_IS_DEFINED_DOUBLE != 0)
  624.     ALIGN_FLOAT (Free);
  625.     Free += 1;
  626. #endif
  627.   Primitive_GC_If_Needed (Length * REAL_SIZE);
  628.   {
  629.     REAL * Temp_Array = ((REAL *) Free);
  630.     Image_Rotate_90cclw (Array, Temp_Array, nrows, ncols);
  631.     C_Array_Copy (Temp_Array, Array, Length);
  632.   }
  633.   {
  634.     SCHEME_OBJECT argument = (ARG_REF (1));
  635.     SET_PAIR_CAR (argument, (LONG_TO_UNSIGNED_FIXNUM (ncols)));
  636.     SET_PAIR_CAR ((PAIR_CDR (argument)), (LONG_TO_UNSIGNED_FIXNUM (nrows)));
  637.     PRIMITIVE_RETURN (argument);
  638.   }
  639. }
  640.  
  641. DEFINE_PRIMITIVE ("IMAGE-MIRROR!", Prim_image_mirror, 1, 1, 0)
  642. {
  643.   long nrows;
  644.   long ncols;
  645.   REAL * Array;
  646.   PRIMITIVE_HEADER (1);
  647.   {
  648.     SCHEME_OBJECT Parray;
  649.     arg_image (1, (&nrows), (&ncols), (&Parray));
  650.     Array = (ARRAY_CONTENTS (Parray));
  651.   }
  652.   C_Mirror_Image (Array, nrows, ncols);    /* side-effecting... */
  653.   PRIMITIVE_RETURN (ARG_REF (1));
  654. }
  655.  
  656.  
  657. /* C routines   referred to above  */
  658.  
  659. /* IMAGE_FAST_TRANSPOSE
  660.    A(i,j) <-> A(j,i) .
  661.    UNWRAP: A(i,j) ----> Array[i*ncols + j]
  662.    convention:= fix row & go by columns .
  663.    UNWRAP is a bijection from the compact plane to the compact interval. */
  664.  
  665. Image_Fast_Transpose (Array, nrows)       /* for square images */
  666.      REAL *Array;
  667.      long nrows;
  668. { long i, j;
  669.   long from, to;
  670.   REAL temp;
  671.   for (i=0;i<nrows;i++) {
  672.     for (j=i;j<nrows;j++) {
  673.       from = i*nrows + j;
  674.       to   = j*nrows + i;    /* (columns transposed-image) = ncols */
  675.       temp        = Array[from];
  676.       Array[from] = Array[to];
  677.       Array[to]   = temp;
  678.     }}
  679. }
  680.  
  681. /* IMAGE_TRANSPOSE
  682.    A(i,j) -> B(j,i) .
  683.    UNWRAP: A(i,j) ----> Array[i*ncols + j]
  684.    convention:= fix row & go by columns .
  685.    UNWRAP is a bijection from the compact plane to the compact interval. */
  686.  
  687. Image_Transpose (Array, New_Array, nrows, ncols)
  688.      REAL *Array, *New_Array;
  689.      long nrows, ncols;
  690. { long i, j;
  691.   for (i=0;i<nrows;i++) {
  692.     for (j=0;j<ncols;j++) {
  693.       /* (columns transposed-image) = nrows */
  694.       New_Array[j*nrows + i] = Array[i*ncols + j];
  695.     }}
  696. }
  697.  
  698. /* IMAGE_ROTATE_90CLW
  699.    A(i,j) <-> A(j, (nrows-1)-i) .
  700.    UNWRAP: A(i,j) ----> Array[i*ncols + j]
  701.    convention:= fix row & go by columns
  702.    UNWRAP is a bijection from the compact plane to the compact interval. */
  703.  
  704. Image_Rotate_90clw (Array, Rotated_Array, nrows, ncols)
  705.      REAL *Array, *Rotated_Array;
  706.      long nrows, ncols;
  707. { long i, j;
  708.  
  709.   for (i=0;i<nrows;i++) {
  710.     for (j=0;j<ncols;j++) {
  711.       /* (columns rotated_image) =nrows */
  712.       Rotated_Array[(j*nrows) + ((nrows-1)-i)] = Array[i*ncols+j];
  713.     }}
  714. }
  715.  
  716. /* ROTATION 90degrees COUNTER-CLOCK-WISE:
  717.    A(i,j) <-> A((nrows-1)-j, i) . (minus 1 because we start from 0).
  718.    UNWRAP: A(i,j) ----> Array[i*ncols + j]
  719.    because of convention:= fix row & go by columns
  720.    UNWRAP is a bijection from the compact plane to the compact interval. */
  721.  
  722. Image_Rotate_90cclw (Array, Rotated_Array, nrows, ncols)
  723.      REAL *Array, *Rotated_Array;
  724.      long nrows, ncols;
  725. { long i, j;
  726.   fast long from_index, to_index;
  727.   long Length=nrows*ncols;
  728.   for (i=0;i<nrows;i++) {
  729.     for (j=0;j<ncols;j++) {
  730.       from_index = i*ncols +j;
  731.       /* (columns rotated-image) = nrows */
  732.       to_index   = ((ncols-1)-j)*nrows + i;
  733.       Rotated_Array[to_index] = Array[from_index];
  734.     }}
  735. }
  736.  
  737. /* IMAGE_MIRROR:
  738.    A(i,j) <-> A(i, (ncols-1)-j)  [ The -1 is there because we count from 0] .
  739.    A(i,j) -------> Array[i*ncols + j]    fix row, read column convention. */
  740.  
  741. C_Mirror_Image (Array, nrows, ncols)
  742.      REAL *Array;
  743.      long nrows, ncols;
  744. { long i, j;
  745.   long ncols2=ncols/2, Length=nrows*ncols;
  746.   REAL temp;
  747.   long from, to;
  748.  
  749.   for (i=0; i<Length; i += ncols) {
  750.     for (j=0; j<ncols2; j++) {    /* DO NOT UNDO the reflections */
  751.       from = i + j;                       /* i is really i*nrows */
  752.       to   = i + (ncols-1)-j;
  753.       temp        = Array[from];
  754.       Array[from] = Array[to];
  755.       Array[to]   = temp;
  756.     }}
  757. }
  758.  
  759. /* IMAGE_ROTATE_90CLW_MIRROR:
  760.    A(i,j) <-> A(j, i)
  761.    this should be identical to image_transpose (see above).
  762.    UNWRAP: A(i,j) ----> Array[i*ncols + j]
  763.    because of convention:= fix row & go by columns
  764.    UNWRAP is a bijection from the compact plane to the compact interval. */
  765.  
  766. C_Rotate_90clw_Mirror_Image (Array, Rotated_Array, nrows, ncols)
  767.      REAL *Array, *Rotated_Array;
  768.      long nrows, ncols;
  769. { long i, j;
  770.   long from, to, Length=nrows*ncols;
  771.  
  772.   for (i=0;i<nrows;i++) {
  773.     for (j=0;j<ncols;j++) {
  774.       from = i*ncols +j;
  775.       /* the columns of the rotated image are nrows! */
  776.       to   = j*nrows +i;
  777.       Rotated_Array[to] = Array[from];
  778.     }}
  779. }
  780.  
  781. DEFINE_PRIMITIVE ("SQUARE-IMAGE-TIME-REVERSE!", Prim_square_image_time_reverse, 2,2, 0)
  782. { long i, rows;
  783.   REAL *a;
  784.   void square_image_time_reverse();
  785.   PRIMITIVE_HEADER (2);
  786.   CHECK_ARG (1, ARRAY_P);
  787.   CHECK_ARG (2, FIXNUM_P);
  788.   a = ARRAY_CONTENTS(ARG_REF(1));
  789.   rows = arg_nonnegative_integer(2);
  790.   if ((rows*rows) != ARRAY_LENGTH(ARG_REF(1)))     error_bad_range_arg(1);
  791.   square_image_time_reverse(a,rows);
  792.  
  793.   PRIMITIVE_RETURN (UNSPECIFIC);
  794. }
  795.  
  796. /* Square Image Time reverse
  797.    is combination of one-dimensional time-reverse
  798.    row-wise and column-wise.
  799.    It can be done slightly more efficiently than below. */
  800.  
  801. void
  802. square_image_time_reverse (x,rows)
  803.      REAL *x;
  804.      long rows;
  805. { long i,cols;
  806.   REAL *xrow, *yrow;
  807.   void C_Array_Time_Reverse();
  808.   cols = rows;            /* square image */
  809.  
  810.   xrow = x;
  811.   for (i=0; i<rows; i++)    /* row-wise */
  812.   { C_Array_Time_Reverse(xrow,cols);
  813.     xrow = xrow + cols; }
  814.  
  815.   Image_Fast_Transpose(x, rows);
  816.  
  817.   xrow = x;
  818.   for (i=0; i<rows; i++)    /* column-wise */
  819.   { C_Array_Time_Reverse(xrow,cols);
  820.     xrow = xrow + cols; }
  821.  
  822.   Image_Fast_Transpose(x, rows);
  823. }
  824.  
  825. /* cs-images */
  826.  
  827. /* operation-1
  828.    groups together procedures     that operate on 1 cs-image-array
  829.    (side-effecting the image) */
  830.  
  831. DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-1!", Prim_cs_image_operation_1, 3,3, 0)
  832. { long rows, opcode;
  833.   REAL *a;
  834.   void cs_image_magnitude(), cs_image_real_part(), cs_image_imag_part();
  835.   PRIMITIVE_HEADER (3);
  836.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  837.   CHECK_ARG (2, FIXNUM_P);    /* rows */
  838.   CHECK_ARG (3, ARRAY_P);    /* input and output image array */
  839.  
  840.   a = ARRAY_CONTENTS(ARG_REF(3));
  841.   rows = arg_nonnegative_integer(2); /*          square images only */
  842.   if ((rows*rows) != ARRAY_LENGTH(ARG_REF(3)))   error_bad_range_arg(1);
  843.   opcode = arg_nonnegative_integer(1);
  844.  
  845.   if (opcode==1)
  846.     cs_image_magnitude(a,rows);
  847.   else if (opcode==2)
  848.     cs_image_real_part(a,rows);
  849.   else if (opcode==3)
  850.     cs_image_imag_part(a,rows);
  851.   else
  852.     error_bad_range_arg(3);    /* illegal opcode */
  853.  
  854.   PRIMITIVE_RETURN (UNSPECIFIC);
  855. }
  856.  
  857.  
  858. void
  859. cs_array_real_part (a,n)
  860.      REAL *a;
  861.      long n;
  862. { long i,n2;            /* works both for even and odd length */
  863.   n2 = n/2;
  864.   for (i=n2+1;i<n;i++) a[i] = a[n-i]; /* copy real values into place */
  865.   /*                                     even signal */
  866. }
  867.  
  868. void
  869. cs_array_imag_part (a,n)
  870.      REAL *a;
  871.      long n;
  872. { long i,n2;
  873.   n2 = n/2;            /* integer division truncates down */
  874.   for (i=n2+1; i<n; i++)    /* works both for even and odd length */
  875.   { a[n-i] = a[i];        /* copy imaginary values into place */
  876.     a[i]   = (-a[i]); }        /* odd signal */
  877.   a[0]     = 0.0;
  878.   if (2*n2 == n)        /* even length, n2 is real only */
  879.     a[n2]    = 0.0;
  880. }
  881.  
  882. /* From now on (below), assume that cs-images (rows=cols)
  883.    have always EVEN LENGTH which is true when they come from FFTs */
  884.  
  885. /*  In the following 3   time-reverse the bottom half rows
  886.     is done to match  the frequencies of complex-images
  887.     coming from cft2d.
  888.     Also transpose is needed to match frequencies identically
  889.  
  890.     #|
  891.     ;; Scrabling of frequencies in  cs-images
  892.  
  893.     ;; start from real image  4x4
  894.  
  895.     ;; rft2d    is a cs-image
  896.     (3.5 .375 -2.75 1.875    -.25 0. 0. -.25    -.25 -.125 0. .125
  897.     .25 .25 0. 0.)
  898.  
  899.     ;; cft2d   transposed
  900.     ;; real
  901.     3.5 .375 -2.75 .375
  902.     -.25  0.  0.  -.25  ; same as cs-image
  903.     -.25 -.125 0. -.125
  904.     -.25 -.25  0.   0.  ; row3 = copy 1 + time-reverse
  905.     ;; imag
  906.     0. 1.875 0. -1.875
  907.     .25 .25 0. 0.       ; same as cs-image
  908.     0. .125 0. -.125
  909.     -.25 0. 0. -.25     ; row 3 = copy 1 + negate + time-reverse
  910.     |#
  911.  
  912.     */
  913.  
  914. void
  915. cs_image_magnitude (x,rows)
  916.      REAL *x;
  917.      long rows;
  918. { long i,j, cols, n,n2, nj; /*     result = real ordinary image */
  919.   REAL *xrow, *yrow;
  920.   cols = rows;            /* input cs-image   is square */
  921.   n = rows;
  922.   n2 = n/2;
  923.  
  924.   xrow = x;
  925.   cs_array_magnitude(xrow, n);  /* row 0 is cs-array */
  926.   xrow = x + n2*cols;
  927.   cs_array_magnitude(xrow, n);  /* row n2 is cs-array */
  928.  
  929.   xrow = x + cols;        /* real part */
  930.   yrow = x + (rows-1)*cols;    /* imag part */
  931.   for (i=1; i<n2; i++) {
  932.     xrow[ 0] = (REAL) sqrt((double) xrow[ 0]*xrow[ 0] + yrow[ 0]*yrow[ 0]);
  933.     xrow[n2] = (REAL) sqrt((double) xrow[n2]*xrow[n2] + yrow[n2]*yrow[n2]);
  934.     yrow[ 0] = xrow[ 0];
  935.     yrow[n2] = xrow[n2];
  936.     for (j=1; j<n2; j++) {
  937.       nj = n-j;
  938.       xrow[ j] = (REAL) sqrt((double) xrow[ j]*xrow[ j] + yrow[ j]*yrow[ j]);
  939.       xrow[nj] = (REAL) sqrt((double) xrow[nj]*xrow[nj] + yrow[nj]*yrow[nj]);
  940.       yrow[j]  = xrow[nj];
  941.       yrow[nj] = xrow[ j];      /* Bottom rows: copy (even) and time-reverse */
  942.     }
  943.     xrow = xrow + cols;
  944.     yrow = yrow - cols; }
  945.   Image_Fast_Transpose(x, n);
  946. }
  947.  
  948. void
  949. cs_image_real_part (x,rows)
  950.      REAL *x;
  951.      long rows;
  952. { long i,j,cols, n,n2;
  953.   REAL *xrow, *yrow;
  954.   void cs_array_real_part();
  955.   cols = rows;            /* square image */
  956.   n = rows;
  957.   n2 = n/2;
  958.  
  959.   xrow = x;
  960.   cs_array_real_part(xrow, n);  /* row 0 is cs-array */
  961.   xrow = x + n2*cols;
  962.   cs_array_real_part(xrow, n);  /* row n2 is cs-array */
  963.  
  964.   xrow = x + cols;        /* real part */
  965.   yrow = x + (rows-1)*cols;    /* imag part */
  966.   for (i=1; i<n2; i++) {
  967.     /* copy real part into imaginary's place (even) */
  968.     yrow[0]  = xrow[0];
  969.     for (j=1; j<n; j++)
  970.       yrow[j] = xrow[n-j];    /* Bottom rows:  copy and time-reverse */
  971.     xrow = xrow + cols;
  972.     yrow = yrow - cols; }
  973.   Image_Fast_Transpose(x, n);
  974. }
  975.  
  976. void
  977. cs_image_imag_part (x,rows)
  978.      REAL *x;
  979.      long rows;
  980. { long i,j,cols, n,n2, nj;
  981.   REAL *xrow, *yrow;
  982.   void cs_array_imag_part();
  983.   cols = rows;            /* square image */
  984.   n = rows;
  985.   n2 = n/2;
  986.  
  987.   xrow = x;
  988.   cs_array_imag_part(xrow, n);  /* row 0 is cs-array */
  989.   xrow = x + n2*cols;
  990.   cs_array_imag_part(xrow, n);  /* row n2 is cs-array */
  991.  
  992.   xrow = x + cols;        /* real part */
  993.   yrow = x + (rows-1)*cols;    /* imag part */
  994.   for (i=1; i<n2; i++) {
  995.     xrow[0]  = yrow[0];        /* copy the imaginary part into real's place */
  996.     xrow[n2] = yrow[n2];
  997.     yrow[0]  = (-yrow[0]);      /* negate (odd) */
  998.     yrow[n2] = (-yrow[n2]);
  999.     for (j=1;j<n2; j++) {
  1000.       nj = n-j;
  1001.       xrow[j]  = yrow[j];     /* copy the imaginary part into real's place */
  1002.       xrow[nj] = yrow[nj];
  1003.       /* Bottom rows: negate (odd) and time-reverse */
  1004.       yrow[j]  = (-xrow[nj]);
  1005.       yrow[nj] = (-xrow[j]); }
  1006.     xrow = xrow + cols;
  1007.     yrow = yrow - cols; }
  1008.   Image_Fast_Transpose(x, n);
  1009. }
  1010.  
  1011. /* cs-image-operation-2
  1012.    groups together procedures     that use 2 cs-image-arrays
  1013.    (usually side-effecting the 2nd image, but not necessarily) */
  1014.  
  1015. DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-2!", Prim_cs_image_operation_2, 4, 4, 0)
  1016. { long rows, nn, opcode;
  1017.   REAL *x,*y;
  1018.   void cs_image_multiply_into_second_one();
  1019.   PRIMITIVE_HEADER (4);
  1020.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1021.   CHECK_ARG (2, FIXNUM_P);    /* rows */
  1022.   CHECK_ARG (3, ARRAY_P);    /* image array 1 */
  1023.   CHECK_ARG (4, ARRAY_P);    /* image array 2 */
  1024.  
  1025.   x = ARRAY_CONTENTS(ARG_REF(3));
  1026.   y = ARRAY_CONTENTS(ARG_REF(4));
  1027.   rows = arg_nonnegative_integer(2); /*          square images only */
  1028.   nn = rows*rows;
  1029.   if (nn != ARRAY_LENGTH(ARG_REF(3)))   error_bad_range_arg(3);
  1030.   if (nn != ARRAY_LENGTH(ARG_REF(4)))   error_bad_range_arg(4);
  1031.  
  1032.   opcode = arg_nonnegative_integer(1);
  1033.  
  1034.   if (opcode==1)
  1035.     cs_image_multiply_into_second_one(x,y,rows); /* result in y */
  1036.   else if (opcode==2)
  1037.     error_bad_range_arg(1);    /* illegal opcode */
  1038.   else
  1039.     error_bad_range_arg(1);    /* illegal opcode */
  1040.  
  1041.   PRIMITIVE_RETURN (UNSPECIFIC);
  1042. }
  1043.  
  1044. void
  1045. cs_image_multiply_into_second_one (x,y, rows)
  1046.      REAL *x,*y;
  1047.      long rows;
  1048. { long i,j,cols, n,n2;
  1049.   REAL *xrow,*yrow,  *xrow_r, *xrow_i, *yrow_r, *yrow_i, temp;
  1050.   cols = rows;            /* square image */
  1051.   n = rows;
  1052.   n2 = n/2;
  1053.  
  1054.   xrow= x; yrow= y;
  1055.   cs_array_multiply_into_second_one(xrow,yrow, n,n2); /*         row 0 */
  1056.  
  1057.   xrow= x+n2*cols; yrow= y+n2*cols;
  1058.   cs_array_multiply_into_second_one(xrow,yrow, n,n2); /*         row n2 */
  1059.  
  1060.   xrow_r= x+cols;           yrow_r= y+cols;
  1061.   xrow_i= x+(n-1)*cols;     yrow_i= y+(n-1)*cols;
  1062.   for (i=1; i<n2; i++) {
  1063.     for (j=0; j<n; j++) {
  1064.       temp      = xrow_r[j]*yrow_r[j]  -  xrow_i[j]*yrow_i[j]; /* real part */
  1065.       yrow_i[j] = xrow_r[j]*yrow_i[j]  +  xrow_i[j]*yrow_r[j]; /* imag part */
  1066.       yrow_r[j] = temp; }
  1067.     xrow_r= xrow_r+cols;   yrow_r= yrow_r+cols;
  1068.     xrow_i= xrow_i-cols;   yrow_i= yrow_i-cols;
  1069.   }
  1070. }
  1071.  
  1072. /* cs-image-operation-2x!     is just like     cs-image-operation-2!
  1073.   but takes an additional flonum argument. */
  1074.  
  1075. DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-2x!", Prim_cs_image_operation_2x, 5, 5, 0)
  1076. { long rows, nn, opcode;
  1077.   REAL *x,*y, flonum_arg;
  1078.   void cs_image_divide_into_z();
  1079.   PRIMITIVE_HEADER (5);
  1080.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1081.   CHECK_ARG (2, FIXNUM_P);    /* rows */
  1082.   CHECK_ARG (3, ARRAY_P);    /* image array 1 */
  1083.   CHECK_ARG (4, ARRAY_P);    /* image array 2 */
  1084.   flonum_arg = (arg_real (5));
  1085.  
  1086.   x = ARRAY_CONTENTS(ARG_REF(3));
  1087.   y = ARRAY_CONTENTS(ARG_REF(4));
  1088.   rows = arg_nonnegative_integer(2); /*          square images only */
  1089.   nn = rows*rows;
  1090.   if (nn != ARRAY_LENGTH(ARG_REF(3)))   error_bad_range_arg(3);
  1091.   if (nn != ARRAY_LENGTH(ARG_REF(4)))   error_bad_range_arg(4);
  1092.  
  1093.   opcode = arg_nonnegative_integer(1);
  1094.  
  1095.   if (opcode==1)
  1096.     cs_image_divide_into_z( x,y, x, rows, flonum_arg); /* result in x */
  1097.   else if (opcode==2)
  1098.     cs_image_divide_into_z( x,y, y, rows, flonum_arg); /* result in y */
  1099.   else
  1100.     error_bad_range_arg(1);    /* illegal opcode */
  1101.  
  1102.   PRIMITIVE_RETURN (UNSPECIFIC);
  1103. }
  1104.  
  1105. /* The convention for inf values in division 1/0 is just like in arrays */
  1106.  
  1107. void
  1108. cs_image_divide_into_z (x,y, z, rows, inf)
  1109.      REAL *x,*y,*z, inf;    /* z can be either x or y */
  1110.      long rows;
  1111. { long i,j,cols, n,n2;
  1112.   REAL temp, radius;
  1113.   REAL  *ar_,*ai_, *br_,*bi_, *zr_,*zi_; /* Letters a,b  correspond to  x,y */
  1114.   REAL *xrow,*yrow,*zrow;
  1115.   cols = rows;            /* square image */
  1116.   n = rows;
  1117.   n2 = n/2;
  1118.  
  1119.   xrow= x; yrow= y; zrow= z;
  1120.   cs_array_divide_into_z( xrow,yrow, zrow, n,n2, inf); /*         row 0 */
  1121.  
  1122.   xrow= x+n2*cols; yrow= y+n2*cols; zrow= z+n2*cols;
  1123.   cs_array_divide_into_z( xrow,yrow, zrow, n,n2, inf); /*         row n2 */
  1124.  
  1125.   ar_= x+cols;           br_= y+cols;            zr_= z+cols;
  1126.   ai_= x+(n-1)*cols;     bi_= y+(n-1)*cols;      zi_= z+(n-1)*cols;
  1127.   for (i=1; i<n2; i++) {
  1128.     for (j=0; j<n; j++) {
  1129.       /* b^2 denominator = real^2 + imag^2 */
  1130.       radius = br_[j]*br_[j]  + bi_[j]*bi_[j];
  1131.  
  1132.       if (radius == 0.0) {
  1133.     if (ar_[j] == 0.0)  zr_[j]  = 1.0;
  1134.     else                zr_[j]  = ar_[j] * inf;
  1135.     if (ai_[j] == 0.0)  zi_[j]  = 1.0;
  1136.     else                zi_[j]  = ai_[j] * inf; }
  1137.       else {
  1138.     temp    =  ar_[j]*br_[j]   +  ai_[j]*bi_[j];
  1139.     zi_[j]  = (ai_[j]*br_[j]   -  ar_[j]*bi_[j]) / radius; /* imag part */
  1140.     zr_[j]  = temp                               / radius; /* real part */
  1141.       }}
  1142.     ar_= ar_+cols;   br_= br_+cols;    zr_= zr_+cols;
  1143.     ai_= ai_-cols;   bi_= bi_-cols;    zi_= zi_-cols;
  1144.   }
  1145. }
  1146.  
  1147. /* operation-3
  1148.    groups together procedures     that use 3 cs-image-arrays
  1149.    (usually side-effecting the 3rd image, but not necessarily) */
  1150.  
  1151. DEFINE_PRIMITIVE ("CS-IMAGE-OPERATION-3!", Prim_cs_image_operation_3, 5, 5, 0)
  1152. { long rows, nn, opcode;
  1153.   REAL *x,*y,*z;
  1154.   void tr_complex_image_to_cs_image();
  1155.   PRIMITIVE_HEADER (5);
  1156.   CHECK_ARG (1, FIXNUM_P);    /* operation opcode */
  1157.   CHECK_ARG (2, FIXNUM_P);    /* rows */
  1158.   CHECK_ARG (3, ARRAY_P);    /* image array 1 */
  1159.   CHECK_ARG (4, ARRAY_P);    /* image array 2 */
  1160.   CHECK_ARG (5, ARRAY_P);    /* image array 3 */
  1161.  
  1162.   x = ARRAY_CONTENTS(ARG_REF(3));
  1163.   y = ARRAY_CONTENTS(ARG_REF(4));
  1164.   z = ARRAY_CONTENTS(ARG_REF(5));
  1165.   rows = arg_nonnegative_integer(2); /*          square images only */
  1166.   nn = rows*rows;
  1167.   if (nn != ARRAY_LENGTH(ARG_REF(3)))   error_bad_range_arg(3);
  1168.   if (nn != ARRAY_LENGTH(ARG_REF(4)))   error_bad_range_arg(4);
  1169.   if (nn != ARRAY_LENGTH(ARG_REF(5)))   error_bad_range_arg(5);
  1170.  
  1171.   opcode = arg_nonnegative_integer(1);
  1172.  
  1173.   if (opcode==1)
  1174.     tr_complex_image_to_cs_image(x,y, z,rows); /* result in z */
  1175.   else if (opcode==2)
  1176.     error_bad_range_arg(1);    /* illegal opcode */
  1177.   else
  1178.     error_bad_range_arg(1);    /* illegal opcode */
  1179.  
  1180.   PRIMITIVE_RETURN (UNSPECIFIC);
  1181. }
  1182.  
  1183. /* x and y must be ALREADY TRANSPOSED real and imaginary parts */
  1184. void
  1185. tr_complex_image_to_cs_image (x,y, z,rows)
  1186.      REAL *x,*y,*z;
  1187.      long rows;
  1188. { long i,j,cols, n,n2, n2_1_n;
  1189.   REAL *xrow, *yrow, *zrow;
  1190.   cols = rows;            /* square image */
  1191.   n = rows;
  1192.   n2 = n/2;
  1193.  
  1194.   xrow= x; yrow= y; zrow= z;
  1195.   for (j=0; j<=n2; j++)
  1196.     /* real part of row 0 (cs-array) */
  1197.     zrow[j] = xrow[j];
  1198.   for (j=n2+1; j<n; j++)
  1199.     /* imag part of row 0 */
  1200.     zrow[j] = yrow[n-j];
  1201.   xrow= x+n2*cols; yrow= y+n2*cols; zrow= z+n2*cols;
  1202.   for (j=0; j<=n2; j++)
  1203.     /* real part of row n2 (cs-array) */
  1204.     zrow[j] = xrow[j];
  1205.   for (j=n2+1; j<n; j++)
  1206.     /* imag part of row n2 */
  1207.     zrow[j] = yrow[n-j];
  1208.   xrow= x+cols;   zrow= z+cols;   n2_1_n = (n2-1)*cols;
  1209.   for (j=0; j<n2_1_n; j++)
  1210.     /* real rows 1,2,..,n2-1 */
  1211.     zrow[j] = xrow[j];
  1212.   yrow= y+(n2-1)*cols;
  1213.   /* imag rows n2+1,n2+2,... */
  1214.   zrow= z+(n2+1)*cols;
  1215.   for (i=1; i<n2; i++) {
  1216.     for (j=0; j<n; j++)   zrow[j] = yrow[j];
  1217.     zrow = zrow + cols;
  1218.     yrow = yrow - cols;
  1219.   }
  1220. }
  1221.