From 5916e93e7ebc112224e28476c3234b6570f35eb5 Mon Sep 17 00:00:00 2001 From: Panayotis Skordos Date: Thu, 22 Jun 1989 21:52:26 +0000 Subject: [PATCH] New 6003 system under construction --- v7/src/microcode/array.c | 1148 ++++++++++++++++++++------------ v7/src/microcode/array.h | 32 +- v7/src/microcode/fft.c | 1357 ++++++++++++++++++++++++++++++++++---- v7/src/microcode/image.c | 1066 +++++++++++++++--------------- 4 files changed, 2522 insertions(+), 1081 deletions(-) diff --git a/v7/src/microcode/array.c b/v7/src/microcode/array.c index 792da02c4..af608f288 100644 --- a/v7/src/microcode/array.c +++ b/v7/src/microcode/array.c @@ -30,20 +30,7 @@ Technology nor of any adaptation thereof in any advertising, promotional, or sales literature without prior written consent from MIT in each case. */ -/* $Header: /Users/cph/tmp/foo/mit-scheme/mit-scheme/v7/src/microcode/Attic/array.c,v 9.35 1989/03/29 04:34:06 pas Exp $ */ - -/* ARRAY = - sequence of REAL(float or double numbers) with a tag on the front */ -/*___________________________________________________________________*/ -/* contents - Scheme_Array datatype - constructors, selectors, operators - procedures for converting between C_Array, and Scheme_Vector - basic and advanced array primitive operations (see manual.scm) */ -/*___________________________________________________________________*/ - -/* See array.h for definition using NM_VECTOR, - and for many useful EXTERN */ +/* $Header: /Users/cph/tmp/foo/mit-scheme/mit-scheme/v7/src/microcode/Attic/array.c,v 9.36 1989/06/22 21:50:18 pas Exp $ */ #include "scheme.h" @@ -53,10 +40,21 @@ MIT in each case. */ #include #include /* contains some math constants */ -/*__________________________________________________________*/ +/* . + ARRAY (as a scheme object) + is a usual array (in C) containing REAL numbers (float/double) + and tagged as a NM_VECTOR + + Basic contents: + constructors, selectors, arithmetic operations, + conversion routines between C_Array, and Scheme_Vector + + see array.h for macros, NM_VECTOR, and extern + */ -/* First some utilities */ + +/* first some utilities */ int Scheme_Number_To_REAL(Arg, Cell) Pointer Arg; REAL *Cell; /* 0 means conversion ok, 1 means too big, 2 means not a number */ @@ -86,7 +84,7 @@ int Scheme_Number_To_REAL(Arg, Cell) Pointer Arg; REAL *Cell; } return (0); } - + int Scheme_Number_To_Double(Arg, Cell) Pointer Arg; double *Cell; /* 0 means conversion ok, 1 means too big, 2 means not a number */ { long Value; @@ -116,15 +114,15 @@ int Scheme_Number_To_Double(Arg, Cell) Pointer Arg; double *Cell; return (0); } -/*__________________begin__________________*/ +/* c */ -/* I think this is not needed, can be done at s-code ... -DEFINE_PRIMITIVE ("ARRAY?", Prim_array_predicate, 1, 1, 0) -{ Primitive_1_Args(); - if (Type_Code(Arg1)==TC_ARRAY) return SHARP_F; - else return SHARP_F; -} -*/ +/* I think this is not needed, it can be done in scheme + DEFINE_PRIMITIVE ("ARRAY?", Prim_array_predicate, 1, 1, 0) + { Primitive_1_Args(); + if (Type_Code(Arg1)==TC_ARRAY) return SHARP_F; + else return SHARP_F; + } + */ DEFINE_PRIMITIVE ("VECTOR->ARRAY", Prim_vector_to_array, 1, 1, 0) { Pointer Scheme_Vector_To_Scheme_Array(); @@ -140,39 +138,41 @@ DEFINE_PRIMITIVE ("ARRAY->VECTOR", Prim_array_to_vector, 1, 1, 0) return Scheme_Array_To_Scheme_Vector(Arg1); } -DEFINE_PRIMITIVE ("ARRAY-CONS", Prim_array_cons, 2, 2, 0) -{ long Length, i, allocated_cells; - REAL Init_Value, *Next; - int Error_Number; +/* array-cons = (array-allocate followed by array-initialize!) + The two are separated because all too often, we only need + array memory space. Even though the initialization is fast, it + happens so often that we get savings. + Also array-initialize! occurs via subarray-offset-scale!. + + */ +DEFINE_PRIMITIVE ("ARRAY-ALLOCATE", Prim_array_allocate, 1,1, 0) +{ long n,allocated_cells; Pointer Result; - Primitive_2_Args(); - Arg_1_Type(TC_FIXNUM); - Range_Check(Length, Arg1, 0, ARRAY_MAX_LENGTH, ERR_ARG_1_BAD_RANGE); + PRIMITIVE_HEADER (1); + CHECK_ARG (1, FIXNUM_P); - Error_Number = Scheme_Number_To_REAL(Arg2, &Init_Value); - if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE); - if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE); + n = arg_nonnegative_integer(1); /* length of array to allocate */ + if (n > ARRAY_MAX_LENGTH) error_bad_range_arg(1); /* avoid memory overflow */ - Allocate_Array(Result,Length,allocated_cells); - Next = Scheme_Array_To_C_Array(Result); - for (i=0; iEnd) Primitive_Error(ERR_ARG_3_BAD_RANGE); - New_Length = (End - Start) + 1; - Allocate_Array(Result, New_Length, allocated_cells); - From_Here = Nth_Array_Loc(Arg1, Start); - To_Here = Scheme_Array_To_C_Array(Result); + +/* ARRAY-COPY! a very powerful primitive + See array.scm for its many applications. + Be Careful when source and destination are the same array. + */ +DEFINE_PRIMITIVE ("SUBARRAY-COPY!", Prim_subarray_copy, 5, 5, 0) +{ long i, i1, i2; + long m, at1, at2; + REAL *a,*b; - C_Array_Copy(From_Here, To_Here, New_Length); - return Result; -} - -DEFINE_PRIMITIVE ("ARRAY-SET-SUBARRAY!", Prim_array_set_subarray, 4, 4, 0) -{ long Length, i, Start, End, New_Length; - REAL *To_Here, *From_Here; - Pointer Result; - - Primitive_4_Args(); - Arg_1_Type(TC_ARRAY); - Arg_2_Type(TC_FIXNUM); - Arg_3_Type(TC_FIXNUM); - Arg_4_Type(TC_ARRAY); - Length = Array_Length(Arg1); - Range_Check(Start, Arg2, 0, Array_Length(Arg1)-1, ERR_ARG_2_BAD_RANGE); - Range_Check(End, Arg3, 0, Array_Length(Arg1)-1, ERR_ARG_3_BAD_RANGE); - if (Start>End) Primitive_Error(ERR_ARG_3_BAD_RANGE); - - New_Length = (End - Start) + 1; - if (New_Length!=Array_Length(Arg4)) Primitive_Error(ERR_ARG_4_BAD_RANGE); - From_Here = Scheme_Array_To_C_Array(Arg4); - To_Here = Nth_Array_Loc(Arg1, Start); + PRIMITIVE_HEADER (5); + CHECK_ARG (1, ARRAY_P); /* source array a */ + CHECK_ARG (2, ARRAY_P); /* destination array b */ - C_Array_Copy(From_Here, To_Here, New_Length); - return Arg1; -} - -DEFINE_PRIMITIVE ("ARRAY-APPEND", Prim_array_append, 2, 2, 0) -{ long Length, Length1, Length2, i, allocated_cells; - REAL *To_Here, *From_Here; - Pointer Result; - - Primitive_2_Args(); - Arg_1_Type(TC_ARRAY); - Arg_2_Type(TC_ARRAY); - Length1 = Array_Length(Arg1); - Length2 = Array_Length(Arg2); - Length = Length1 + Length2; - - Allocate_Array(Result, Length, allocated_cells); - To_Here = Scheme_Array_To_C_Array(Result); - From_Here = Scheme_Array_To_C_Array(Arg1); - - for (i=0; i < Length1; i++) { - *To_Here++ = *From_Here; - From_Here++ ; - } + a = Scheme_Array_To_C_Array(ARG_REF(1)); + b = Scheme_Array_To_C_Array(ARG_REF(2)); + at1 = arg_nonnegative_integer(3); /* at1 = starting index in source array */ + at2 = arg_nonnegative_integer(4); /* at2 = starting index in destination array */ + m = arg_nonnegative_integer(5); /* m = number of points to copy */ + + if ((at1 + m) > (Array_Length(ARG_REF(1)))) error_bad_range_arg(3); + if ((at2 + m) > (Array_Length(ARG_REF(2)))) error_bad_range_arg(4); + /* These 2 checks cover all cases */ - From_Here = Scheme_Array_To_C_Array(Arg2); - for (i=0; i < Length2; i++) { - *To_Here++ = *From_Here; - From_Here++ ; - } + for (i=0,i1=at1,i2=at2; i (Array_Length(ARG_REF(1)))) error_bad_range_arg(3); - Result = Arg1; - From_Here = Scheme_Array_To_C_Array(Arg1); - To_Here = Scheme_Array_To_C_Array(Result); + errcode = Scheme_Number_To_REAL(ARG_REF(4), &offset); + if (errcode==1) error_bad_range_arg(4); if (errcode==2) error_wrong_type_arg(4); + errcode = Scheme_Number_To_REAL(ARG_REF(5), &scale); + if (errcode==1) error_bad_range_arg(5); if (errcode==2) error_wrong_type_arg(5); + + if ((offset == 0.0) && (scale == 1.0)) + ; /* be smart */ + else if (scale == 0.0) + for (i=at; i (Array_Length(ARG_REF(1)))) error_bad_range_arg(4); + + errcode = Scheme_Number_To_Double(ARG_REF(5), &x); + if (errcode==1) error_bad_range_arg(5); if (errcode==2) error_wrong_type_arg(5); + errcode = Scheme_Number_To_Double(ARG_REF(6), &y); + if (errcode==1) error_bad_range_arg(6); if (errcode==2) error_wrong_type_arg(6); + + a = Scheme_Array_To_C_Array(ARG_REF(1)); + b = Scheme_Array_To_C_Array(ARG_REF(2)); + if ((Array_Length(ARG_REF(1))) != (Array_Length(ARG_REF(2)))) error_bad_range_arg(2); + + if (x==0.0) /* imaginary only */ + if (y==0.0) + for (i=at; i= 0.0) /* It may be faster to look at the sign of mantissa and dispatch */ + if ((*a) >= 0.0) /* It may be faster to look at the sign + of mantissa, and dispatch */ modf( ((double) ((*a)+0.5)), &integral_part); else modf( ((double) ((*a)-0.5)), &integral_part); @@ -541,8 +757,9 @@ void REALbessel2(order,a,b) long order; REAL *a,*b; /* Bessel of second kind */ (*b) = ( (REAL) yn(((int) order), ((double) (*a))) ); } -/* Table to store the available functions for transforming arrays. - It also stores the corresponding numofargs (whether unary or binary function). +/* Table to store the available unary-functions. + Also some binary functions at the end -- not available right now. + The (1 and 2)s denote the numofargs (1 for unary 2 for binary) */ struct array_func_table { @@ -597,32 +814,39 @@ DEFINE_PRIMITIVE ("ARRAY-UNARY-FUNCTION!", Prim_array_unary_function, 2, 2, 0) return Result; } -/* The following is accumulate of + and * - code numbers are 0 1 + +/* Accumulate + using combinators + or * + corresponding type codes 0 1 */ -DEFINE_PRIMITIVE ("ARRAY-ACCUMULATE", Prim_array_accumulate, 2, 2, 0) -{ long Length, i; - REAL *a, result; - long functc; +DEFINE_PRIMITIVE ("SUBARRAY-ACCUMULATE", Prim_subarray_accumulate, 4,4, 0) +{ long at,m,mplus, tc, i; + REAL *a; + double result; - Primitive_2_Args(); - Arg_1_Type(TC_ARRAY); - Arg_2_Type(TC_FIXNUM); - Length = Array_Length(Arg1); - Range_Check(functc, Arg2, 0, 1, ERR_ARG_2_BAD_RANGE); - a = Scheme_Array_To_C_Array(Arg1); + PRIMITIVE_HEADER (4); + CHECK_ARG (1, ARRAY_P); /* a = input array */ + a = Scheme_Array_To_C_Array(ARG_REF(1)); + tc = arg_nonnegative_integer(2); /* tc = type code 0 or 1 */ + at = arg_nonnegative_integer(3); /* at = starting index */ + m = arg_nonnegative_integer(4); /* m = number of points to process */ - if (functc==0) + mplus = at + m; + if (mplus > (Array_Length(ARG_REF(1)))) error_bad_range_arg(4); + + if (tc==0) { result = 0.0; - for (i=0;i (Array_Length(ARG_REF(1)))) error_bad_range_arg(3); + + C_Array_Find_Min_Max ( &(a[at]), m, &nmin, &nmax); + nmin = nmin + at; /* offset appropriately */ + nmax = nmax + at; + Primitive_GC_If_Needed(4); Result = Make_Pointer(TC_LIST, Free); Orig_Free = Free; @@ -673,8 +906,10 @@ DEFINE_PRIMITIVE ("ARRAY-MIN-MAX-INDEX", Prim_array_min_max_index, 1, 1, 0) *Orig_Free++ = Make_Pointer(TC_LIST, Orig_Free+1); *Orig_Free++ = Make_Non_Pointer(TC_FIXNUM, nmax); *Orig_Free=EMPTY_LIST; - return Result; + + PRIMITIVE_RETURN (Result); } + void C_Array_Find_Min_Max(x, n, nmin, nmax) REAL *x; long n, *nmax, *nmin; { REAL *xold = x; register REAL xmin, xmax; @@ -705,11 +940,14 @@ void C_Array_Find_Min_Max(x, n, nmin, nmax) REAL *x; long n, *nmax, *nmin; } -/* The following becomes obsolete. - Done using array-reduce + divide by array-length +/* array-average + can be done with (array-reduce +) and division by array-length. + But there is also this C primitive. Keep it around, may be useful someday. */ + DEFINE_PRIMITIVE ("ARRAY-AVERAGE", Prim_array_find_average, 1, 1, 0) { long Length; REAL average; + void C_Array_Find_Average(); Primitive_1_Args(); Arg_1_Type(TC_ARRAY); Length = Array_Length(Arg1); @@ -717,10 +955,12 @@ DEFINE_PRIMITIVE ("ARRAY-AVERAGE", Prim_array_find_average, 1, 1, 0) C_Array_Find_Average( Scheme_Array_To_C_Array(Arg1), Length, &average); Reduced_Flonum_Result((double) average); } + /* Computes the average in pieces, so as to reduce roundoff smearing in cumulative sum. example= first huge positive numbers, then small nums, then huge negative numbers. */ + void C_Array_Find_Average(Array, Length, pAverage) long Length; REAL *Array, *pAverage; { long i; @@ -772,23 +1012,22 @@ void C_Array_Make_Histogram(Array, Length, Histogram, npoints) } DEFINE_PRIMITIVE ("ARRAY-CLIP-MIN-MAX!", Prim_array_clip_min_max, 3, 3, 0) -{ long Length, i; /* , allocated_cells; */ +{ long Length, i; REAL *To_Here, *From_Here, xmin, xmax; - Pointer Result; - int Error_Number; + int errcode; Primitive_3_Args(); Arg_1_Type(TC_ARRAY); - Error_Number=Scheme_Number_To_REAL(Arg2, &xmin); - if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE); - if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE); - Error_Number=Scheme_Number_To_REAL(Arg3, &xmax); - if (Error_Number == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE); - if (Error_Number == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE); + errcode=Scheme_Number_To_REAL(Arg2, &xmin); + if (errcode == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE); + if (errcode == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE); + errcode=Scheme_Number_To_REAL(Arg3, &xmax); + if (errcode == 1) Primitive_Error(ERR_ARG_3_BAD_RANGE); + if (errcode == 2) Primitive_Error(ERR_ARG_3_WRONG_TYPE); Length = Array_Length(Arg1); - Result = Arg1; + From_Here = Scheme_Array_To_C_Array(Arg1); - To_Here = Scheme_Array_To_C_Array(Result); + To_Here = Scheme_Array_To_C_Array(Arg1); if (xmin>xmax) Primitive_Error(ERR_ARG_3_BAD_RANGE); for (i=0; i < Length; i++) { @@ -796,15 +1035,16 @@ DEFINE_PRIMITIVE ("ARRAY-CLIP-MIN-MAX!", Prim_array_clip_min_max, 3, 3, 0) else if ((*From_Here)>xmax) *To_Here++ = xmax; else *To_Here++ = *From_Here; From_Here++ ; } - return Result; + return SHARP_F; } -DEFINE_PRIMITIVE ("ARRAY-MAKE-POLAR!", Prim_array_make_polar, 2, 2, 0) +DEFINE_PRIMITIVE ("COMPLEX-ARRAY-TO-POLAR!", + Prim_complex_array_to_polar, 2, 2, 0) { long Length, i; REAL *To_Here_Mag, *To_Here_Phase; REAL *From_Here_Real, *From_Here_Imag; - Pointer Result_Mag, Result_Phase, answer; - + Pointer Result_Mag, Result_Phase; + Primitive_2_Args(); Arg_1_Type(TC_ARRAY); Arg_2_Type(TC_ARRAY); @@ -820,48 +1060,163 @@ DEFINE_PRIMITIVE ("ARRAY-MAKE-POLAR!", Prim_array_make_polar, 2, 2, 0) To_Here_Phase = Scheme_Array_To_C_Array(Result_Phase); for (i=0; i < Length; i++) - { - C_Make_Polar(*From_Here_Real, *From_Here_Imag, *To_Here_Mag, *To_Here_Phase); + { C_Make_Polar(*From_Here_Real, *From_Here_Imag, *To_Here_Mag, *To_Here_Phase); From_Here_Real++ ; From_Here_Imag++ ; - To_Here_Mag++ ; To_Here_Phase++ ; - } + To_Here_Mag++ ; To_Here_Phase++ ; } - Primitive_GC_If_Needed(4); - answer = Make_Pointer(TC_LIST, Free); - *Free++ = Result_Mag; - *Free = Make_Pointer(TC_LIST, Free+1); - Free += 1; - *Free++ = Result_Phase; - *Free++ = EMPTY_LIST; - return answer; + return SHARP_F; } -DEFINE_PRIMITIVE ("ARRAY-FIND-MAGNITUDE", Prim_array_find_magnitude, 2, 2, 0) -{ long Length, i, allocated_cells; - REAL *From_Here_Real, *From_Here_Imag, *To_Here; - Pointer Result; +DEFINE_PRIMITIVE ("COMPLEX-ARRAY-MAGNITUDE!", Prim_complex_array_magnitude, 3, 3, 0) +{ long n, i; + REAL *a, *b, *mg; - Primitive_2_Args(); - Arg_1_Type(TC_ARRAY); - Arg_2_Type(TC_ARRAY); - Length = Array_Length(Arg1); - if (Length != Array_Length(Arg2)) Primitive_Error(ERR_ARG_1_BAD_RANGE); + PRIMITIVE_HEADER (3); + CHECK_ARG (1, ARRAY_P); /* input array -- n real part */ + CHECK_ARG (2, ARRAY_P); /* input array -- n imag part */ + CHECK_ARG (3, ARRAY_P); /* ouput array -- n magnitude */ + + n = Array_Length(ARG_REF(1)); + if (n != Array_Length(ARG_REF(2))) error_bad_range_arg(2); + if (n != Array_Length(ARG_REF(3))) error_bad_range_arg(3); + + a = Scheme_Array_To_C_Array(ARG_REF(1)); /* real part */ + b = Scheme_Array_To_C_Array(ARG_REF(2)); /* imag part */ + mg = Scheme_Array_To_C_Array(ARG_REF(3)); /* magnitude */ + + for (i=0; i information is lost when going from rect to polar (compromise to save space) + + To invert from polar to rect, assume the signal is a continuous curve + and choose sign (from angle 0 or pi) at index i=0 same sign as i=1, + and i=n2 same sign as i=n2+1 + */ + +DEFINE_PRIMITIVE ("CS-ARRAY-TO-POLAR!", Prim_cs_array_to_polar, 1,1, 0) +{ long n, i; + REAL *a; + void cs_array_to_polar(); + + PRIMITIVE_HEADER (1); + CHECK_ARG (1, ARRAY_P); /* input and output array -- both cs-arrays */ + a = Scheme_Array_To_C_Array(ARG_REF(1)); + n = Array_Length(ARG_REF(1)); + + cs_array_to_polar(a,n); + PRIMITIVE_RETURN (NIL); +} + +void cs_array_to_polar(a,n) + REAL *a; long n; +{ long i, n2; + double real, imag; /* temporary variables */ + n2 = n/2; /* integer division truncates down */ + + a[0] = (REAL) fabs((double) a[0]); /* implicitly angle = 0, but it could be pi */ + if (2*n2 == n) /* even length, n2 is only real */ + a[n2] = (REAL) fabs((double) a[n2]); /* implicitly angle = 0, but it could be pi */ + else + n2 = n2+1; /* odd length, make the loop include the n2 index */ + + for (i=1; i 0.0) + a[0] = a[0]; /* assume angle = 0 */ + else + a[0] = (-a[0]); /* assume angle = pi */ + + if (2*n2 == n) /* even length, n2 is real only */ + if (a[n2+1] > 0.0) + a[n2] = a[n2]; /* assume angle = 0 */ + else + a[n2] = (-a[n2]); /* assume angle = pi */ + else + n2 = n2+1; /* odd length, make the loop include the n2 index */ + + for (i=1; i1; power++) + { if ( (i % 2) == 1) error_bad_range_arg(2); + i=i/2; } + + f1 = Scheme_Array_To_C_Array(ARG_REF(2)); + f2 = Scheme_Array_To_C_Array(ARG_REF(3)); + if (f1==f2) error_wrong_type_arg(2); + + wcos = Scheme_Array_To_C_Array(ARG_REF(4)); /* twiddle tables */ + if (Array_Length(ARG_REF(4)) != (3*length/4)) error_bad_range_arg(4); + w3cos = wcos + length/4; + w3sin = w3cos + length/4; + if ((arg_nonnegative_integer(5)) == 1) + pas_cft(1, flag, f1,f2, length, power, wcos,w3cos,w3sin); + else + pas_cft(0, flag, f1,f2, length, power, wcos,w3cos,w3sin); + /* 1 means tables are already made + 0 means compute new tables */ + + + PRIMITIVE_RETURN (NIL); +} + +DEFINE_PRIMITIVE ("PAS-CFT-MAKE-TWIDDLE-TABLES!", + Prim_pas_cft_make_twiddle_tables, 2, 2, 0) +{ long length, power, i; + REAL *wcos,*w3cos,*w3sin; + void pas_cft_make_twiddle_tables_once(); + PRIMITIVE_HEADER (2); + + length = arg_nonnegative_integer(1); /* length of cft that we intend to compute */ + CHECK_ARG (2, ARRAY_P); /* storage for twiddle tables */ + if (Array_Length(ARG_REF(2)) != (3*length/4)) error_bad_range_arg(2); + + power=0; + for (power=0, i=length; i>1; power++) + { if ( (i % 2) == 1) error_bad_range_arg(1); + i=i/2; } + + wcos = Scheme_Array_To_C_Array(ARG_REF(2)); /* twiddle tables */ + w3cos = wcos + length/4; + w3sin = w3cos + length/4; + pas_cft_make_twiddle_tables_once(length, power, wcos,w3cos,w3sin); + + PRIMITIVE_RETURN (NIL); +} + + +/* + C COMPLEX FOURIER TRANSFORM (Split-Radix, Decimation-in-frequency) + C (adapted and optimized from Sorensen,et.al. ASSP-34 no.1 page 152, February 1986) + */ + +/* Twiddle Tables for PAS_CFT; + (tables for forward transform only) + Inverse transform === forward CFT (without 1/N scaling) followed by time-reversal. + / + The tables contain (2pi/N)*i for i=0,1,2,..,N/4 + (except i=0 is ignored, never used) + / + Table for wsin[i] is not needed because wsin[i]=wcos[n4-i]. + Table for w3sin[i] is needed however. The previous relationship does not work for w3sin. + */ + +/* There are two routines for making twiddle tables: + a fast one, and a slower one but more precise. + The differences in speed and accuracy are actually rather small, but anyway. + Use the slow one for making permanent tables. + */ + +void pas_cft_make_twiddle_tables(n,m, wcos,w3cos,w3sin) /* efficient version */ + REAL *wcos, *w3cos, *w3sin; + long n,m; +{ long i, n4; + double tm; + REAL costm,sintm; + n4 = n/4; + for (i=1; i= 4 */ + REAL *x,*y, *wcos,*w3cos,*w3sin; + long n,m; +{ /* REAL a,a3,e; no need anymore, use tables */ + REAL r1,r2,s1,s2,s3, xt, cc1,cc3,ss1,ss3; + long n1,n2,n4, i,j,k, is,id, i0,i1,i2,i3; + long windex0, windex, windex_n4; /* indices for twiddle tables */ + /********** fortran indices start from 1,... **/ + x = x-1; /* TRICK---- x(0) is now illegal, but x(1) and x(n) are valid */ + y = y-1; + /********** fortran indices start from 1,... **/ + /* c */ + /* c-----first M-1 stages of transform */ + /* c */ + windex_n4 = n/4; /* need for indexing sin via wcos twiddle table */ + n2 = 2*n; + for (k=1; k>1; /* n2 = n2/2; */ + n4 = n2>>2; /* n4 = n2/4; */ + /* e = 6.283185307179586476925287 / ((REAL) n2); no need anymore, use tables */ + /* a = 0.0; */ + { + /* j=1; */ + /* The first iteration in the loop "DO 20 J = 1, N4" + is done specially to save operations involving sin=0, cos=1 */ + /* a = j*e; no need anymore, use tables */ + is = 1; /* is = j; */ + id = 2*n2; + label40first: + for (i0=is; i0= j) goto label101; /* if (i .ge. j) goto 101 */ + xt = x[j]; + x[j] = x[i]; + x[i] = xt; + xt = y[j]; + y[j] = y[i]; + y[i] = xt; + label101: k = n>>1; /* k = n/2; */ + label102: if (k>=j) goto label103; + j = j - k; + k = k>>1; /* k = k/2; */ + goto label102; + label103: j = j + k; + } /* 104 CONTINUE */ + /* c-------------------------------------*/ + /* c */ +} /* RETURN END */ + + + +DEFINE_PRIMITIVE ("PAS-RFT-CSFT!", Prim_pas_rft_csft, 5, 5, 0) +{ long i, length, power, flag, ft_type; + REAL *f1, *wcos,*w3cos,*w3sin; + void pas_rft(), pas_csft(); + PRIMITIVE_HEADER (5); + CHECK_ARG (1, FIXNUM_P); /* flag 1=forward, else backward transform */ + CHECK_ARG (2, ARRAY_P); /* Input data (real or cs) */ + CHECK_ARG (3, ARRAY_P); /* Twiddle tables, total length = 4*(length/8) */ + CHECK_ARG (4, FIXNUM_P); /* (1)=tables precomputed, else recompute */ + CHECK_ARG (5, FIXNUM_P); /* ft_type = 1 or 3 + 1 means compute rft, 3 means compute csft */ + flag = Get_Integer(ARG_REF(1)); + f1 = Scheme_Array_To_C_Array(ARG_REF(2)); + length = Array_Length(ARG_REF(2)); + for (power=0, i=length; i>1; power++) + { if ( (i % 2) == 1) error_bad_range_arg(2); + i=i/2; } + + wcos = Scheme_Array_To_C_Array(ARG_REF(3)); /* twiddle tables */ + if (Array_Length(ARG_REF(3)) != (4*length/8)) error_bad_range_arg(3); + w3cos = wcos + (length/4); + w3sin = w3cos + (length/8); + + ft_type = (arg_nonnegative_integer(5)); /* rft or csft */ + if (ft_type == 1) { + if ((arg_nonnegative_integer(4)) == 1) + pas_rft (1, flag, f1, length, power, wcos,w3cos,w3sin); + else pas_rft (0, flag, f1, length, power, wcos,w3cos,w3sin); + } + else if (ft_type == 3) { + if ((arg_nonnegative_integer(4)) == 1) + pas_csft (1, flag, f1, length, power, wcos,w3cos,w3sin); + else pas_csft (0, flag, f1, length, power, wcos,w3cos,w3sin); + /* 1 means tables are already made + 0 means compute new tables */ + } + else error_bad_range_arg(5); + + PRIMITIVE_RETURN (NIL); +} + +DEFINE_PRIMITIVE ("PAS-REALDATA-MAKE-TWIDDLE-TABLES!", + Prim_pas_realdata_make_twiddle_tables, 2, 2, 0) +{ long length, power, i; + REAL *wcos,*w3cos,*w3sin; + void pas_realdata_make_twiddle_tables_once(); + PRIMITIVE_HEADER (2); + + length = arg_nonnegative_integer(1); /* length of rft that we intend to compute */ + CHECK_ARG (2, ARRAY_P); /* storage for twiddle tables */ + if (Array_Length(ARG_REF(2)) != (4*length/8)) error_bad_range_arg(2); + + power=0; + for (power=0, i=length; i>1; power++) + { if ( (i % 2) == 1) error_bad_range_arg(1); + i=i/2; } + + wcos = Scheme_Array_To_C_Array(ARG_REF(2)); /* twiddle tables */ + w3cos = wcos + length/4; + w3sin = w3cos + length/8; + pas_realdata_make_twiddle_tables_once(length, power, wcos,w3cos,w3sin); + + PRIMITIVE_RETURN (NIL); +} + +/* + C REAL FOURIER TRANSFORM (Split-Radix, Decimation-in-time) + C (adapted from Sorensen,et.al. ASSP-35 no.6 page 849, October 1986) + C + C the output is [Re(0),Re(1),...,Re(n/2), Im(n/2-1),...,Im(1)] + */ + +/* Twiddle Tables for PAS_RFT and PAS_CSFT + are identical. -> pas_realdata_make_twiddle_tables + (but they are indexed differently in each case) + / + The tables contain (2pi/N)*i where i=0,1,2,..,N/4 for wcos + and (2pi/N)*i where i=0,1,2,..,N/8 for w3cos w3sin + (except i=0 is ignored, never used) + / + Table for wsin[i] is not needed because wsin[i]=wcos[n4-i]. + Table for w3sin[i] is needed however. The previous relationship does not work for w3sin. + / + Instead of getting sin() from a wsin[i] table with i=1,..,N/8 + we get it from wcos[n4-i]. + This way we can use a CFT table which goes up to N/4 + for RFT CSFT also. We do so in image-processing (rft2d-csft2d). + */ + +/* There are two routines for making twiddle tables: + a fast one, and a slower one but more precise. + The differences in speed and accuracy are actually rather small, but anyway. + Use the slow one for making tables that stay around. + */ + +void pas_realdata_make_twiddle_tables(n,m, wcos,w3cos,w3sin) /* efficient version */ + REAL *wcos, *w3cos, *w3sin; + long n,m; +{ long i, n4, n8; + double tm; + REAL costm,sintm; + n4 = n/4; + n8 = n/8; + for (i=1; i= j) goto label101; /* if (i .ge. j) goto 101 */ + xt = x[j]; + x[j] = x[i]; + x[i] = xt; + label101: k = n>>1; /* k = n/2; */ + label102: if (k>=j) goto label103; + j = j - k; + k = k>>1; /* k = k/2; */ + goto label102; + label103: j = j + k; + } /* 104 CONTINUE */ + /* c-------------------------------------*/ + /* c */ + /* c ----length-two-butterflies----------- */ + is = 1; + id = 4; + label70: + for (i0=is; i0<=n; i0=i0+id) /* 70 DO 60 I0 = IS,N,ID */ + { i1 = i0 + 1; + r1 = x[i0]; + x[i0] = r1 + x[i1]; + x[i1] = r1 - x[i1]; + } /* 60 CONTINUE */ + is = 2*id - 1; + id = 4*id; + if (is < n) goto label70; /* IF (IS.LT.N) GOTO 70 */ + /* C + C -------L-shaped-butterflies-------- */ + n2 = 2; + for (k=2; k<=m; k++) /* DO 10 K = 2,M */ + { n2 = n2 * 2; + n4 = n2>>2; /* n4 = n2/4; */ + n8 = n2>>3; /* n8 = n2/8; */ + /* e = 6.283185307179586476925287 / ((REAL) n2); no need anymore, use tables */ + is = 0; + id = n2 * 2; + label40: + for (i=is; i>1; /* n2 = n2/2; */ + n4 = n2>>2; /* n4 = n2/4; */ + n8 = n4>>1; /* n8 = n4/2; */ + /* e = 6.283185307179586476925287 / ((REAL) n2); no need anymore, use tables */ + label17: + for (i=is; i= j) goto label101; /* if (i .ge. j) goto 101 */ + xt = x[j]; + x[j] = x[i]; + x[i] = xt; + label101: k = n>>1; /* k = n/2; */ + label102: if (k>=j) goto label103; + j = j - k; + k = k>>1; /* k = k/2; */ + goto label102; + label103: j = j + k; + } /* 104 CONTINUE */ + /* c */ +} /* RETURN END */ + + + + +/* Image processing only for square images + (old stuff handles non-square but is slow) + For 2d FTs precomputed tables or not, make almost no difference in total time. + */ + +DEFINE_PRIMITIVE ("PAS-CFT2D!", Prim_pas_cft2d, 5,5, 0) +{ long i, length, power, flag, rows,rowpower; + REAL *f1,*f2, *wcos,*w3cos,*w3sin; + void pas_cft2d(); + PRIMITIVE_HEADER (5); + CHECK_ARG (1, FIXNUM_P); /* flag forward-backward */ + CHECK_ARG (2, ARRAY_P); /* real part */ + CHECK_ARG (3, ARRAY_P); /* imag part */ + CHECK_ARG (4, ARRAY_P); /* twiddle tables, length = 3*(rows/4) */ + CHECK_ARG (5, FIXNUM_P); /* (1)=tables precomputed, else recompute */ + + flag = Get_Integer(ARG_REF(1)); + length = Array_Length(ARG_REF(2)); + if (length != (Array_Length(ARG_REF(3)))) error_bad_range_arg(2); + + for (power=0, i=length; i>1; power++) /* length must be power of 2 */ + { if ( (i % 2) == 1) error_bad_range_arg(2); + i=i/2; } + + if ((power % 2) == 1) error_bad_range_arg(2); + rowpower = (power/2); + rows = (1<1; power++) /* length must be power of 2 */ + { if ( (i % 2) == 1) error_bad_range_arg(2); + i=i/2; } + + if ((power % 2) == 1) error_bad_range_arg(2); + rowpower = (power/2); + rows = (1<>1, a; @@ -176,33 +1303,13 @@ void C_Array_FFT_With_Given_Tables(flag, power, n, f1, f2, g1,g2,w1,w2) } -/* CHIRP-Z-TRANSFORM (for complex data) +/* CHIRP-Z-TRANSFORM (complex data) */ -#define take_modulo_one(x, answer) \ -{ long ignore_integral_part; \ - double modf(); \ - answer = (REAL) modf( ((double) x), &ignore_integral_part); } - -#define power_of_2_to_power(n,power) \ -{ long i; \ - for (power=0,i=n; i>1; power++) \ - { if ((i%2) == 1) { printf("\nNot a power of 2--- %d\n",n); exit(1); } \ - i=i/2; }} - -long smallest_power_of_2_ge(n) - long n; -{ long i,power; - if (n<0) { printf("\nsmallest_pwr_of_2_ge negative argument--- %d\n", n); exit(1); } - power=0; i=1; - while (i1; power++) { - if ( (i % 2) == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE); + if ( (i % 2) == 1) error_bad_range_arg(2); i=i/2; } - f1 = Scheme_Array_To_C_Array(Arg2); - f2 = Scheme_Array_To_C_Array(Arg3); - if (f1==f2) Primitive_Error(ERR_ARG_2_WRONG_TYPE); + f1 = Scheme_Array_To_C_Array(ARG_REF(2)); + f2 = Scheme_Array_To_C_Array(ARG_REF(3)); + if (f1==f2) error_wrong_type_arg(2); Primitive_GC_If_Needed(length*3*REAL_SIZE); Work_Here = (REAL *) Free; @@ -492,75 +1615,67 @@ DEFINE_PRIMITIVE ("ARRAY-FFT!", Prim_array_fft, 3, 3, 0) C_Array_FFT(flag, power, length, f1,f2,g1,g2,w1,w2); - Primitive_GC_If_Needed(4); - answer = Make_Pointer(TC_LIST, Free); - *Free++ = Arg2; - *Free = Make_Pointer(TC_LIST, Free+1); - Free += 1; - *Free++ = Arg3; - *Free++ = NIL; - return answer; + PRIMITIVE_RETURN (NIL); } -DEFINE_PRIMITIVE ("ARRAY-CZT", Prim_array_czt, 5, 5, 0) +DEFINE_PRIMITIVE ("ARRAY-CZT!", Prim_array_czt, 6,6, 0) { double phi,rho; - long N,M,L; - long log2_L,maxMN,smallest_power_of_2_ge(), allocated_cells; - int Error_Number, Scheme_Number_To_Double(); - Pointer answer, answer_1,answer_2; - REAL *f1,*f2,*fo1,*fo2, *g1,*g2, *fft_w1,*fft_w2,*czt_w1,*czt_w2; - REAL *Work_Here; - Primitive_5_Args(); - Arg_3_Type(TC_FIXNUM); - Arg_4_Type(TC_ARRAY); /* real */ - Arg_5_Type(TC_ARRAY); /* imag */ - Set_Time_Zone(Zone_Math); - - Error_Number = Scheme_Number_To_Double(Arg1, &phi); /* phi=starting point [0,1]*/ - if (Error_Number == 1) Primitive_Error(ERR_ARG_1_BAD_RANGE); - if (Error_Number == 2) Primitive_Error(ERR_ARG_1_WRONG_TYPE); - Error_Number = Scheme_Number_To_Double(Arg2, &rho); /* rho=resolution [0,1] */ - if (Error_Number == 1) Primitive_Error(ERR_ARG_2_BAD_RANGE); - if (Error_Number == 2) Primitive_Error(ERR_ARG_2_WRONG_TYPE); - - M = Get_Integer(Arg3); /* M=length of output data */ - N = Array_Length(Arg4); /* N=length of input data */ - if (N != (Array_Length(Arg5))) Primitive_Error(ERR_ARG_5_BAD_RANGE); - if ((M+N-1)<4) Primitive_Error(ERR_ARG_3_BAD_RANGE); + long N,M,L, i; + long log2_L, maxMN; + long smallest_power_of_2_ge(); + int errcode; + REAL *a,*b,*c,*d; + REAL *f1,*f2,*fo1,*fo2, *g1,*g2, *fft_w1,*fft_w2,*czt_w1,*czt_w2, *Work_Here; + + PRIMITIVE_HEADER (6); + /* 1 phi=starting point [0,1]*/ + /* 2 rho=resolution [0,1] */ + CHECK_ARG (3, ARRAY_P); /* input real */ + CHECK_ARG (4, ARRAY_P); /* input imag */ + CHECK_ARG (5, ARRAY_P); /* output real */ + CHECK_ARG (6, ARRAY_P); /* output imag */ + + errcode = Scheme_Number_To_Double(ARG_REF(1), &phi); + if (errcode==1) error_bad_range_arg(1); if (errcode==2) error_wrong_type_arg(1); + errcode = Scheme_Number_To_Double(ARG_REF(2), &rho); + if (errcode==1) error_bad_range_arg(2); if (errcode==2) error_wrong_type_arg(2); + + a = Scheme_Array_To_C_Array(ARG_REF(3)); + b = Scheme_Array_To_C_Array(ARG_REF(4)); + c = Scheme_Array_To_C_Array(ARG_REF(5)); + d = Scheme_Array_To_C_Array(ARG_REF(6)); + + N = Array_Length(ARG_REF(3)); /* N = input length */ + M = Array_Length(ARG_REF(5)); /* M = output length */ + if (N!=(Array_Length(ARG_REF(4)))) error_bad_range_arg(3); + if (M!=(Array_Length(ARG_REF(6)))) error_bad_range_arg(5); + + if ((M+N-1) < 4) error_bad_range_arg(5); log2_L = smallest_power_of_2_ge(M+N-1); L = 1<Array[i*ncols+j] */ - }} - - return Result; -} - -/* The following does not work properly, to be fixed if need. +/* The following does not work, to be fixed. */ DEFINE_PRIMITIVE ("IMAGE-DOUBLE-TO-FLOAT!", Prim_image_double_to_float, 1, 1, 0) { long Length; @@ -434,7 +377,6 @@ DEFINE_PRIMITIVE ("IMAGE-DOUBLE-TO-FLOAT!", Prim_image_double_to_float, 1, 1, 0) *To_Here = ((float) temp_value_cell); To_Here++; } - /* and now SIDE-EFFECT the ARRAY_HEADER */ allocated_cells = (Length * ((sizeof(Pointer)+sizeof(float)-1) / sizeof(Pointer)) + @@ -446,157 +388,137 @@ DEFINE_PRIMITIVE ("IMAGE-DOUBLE-TO-FLOAT!", Prim_image_double_to_float, 1, 1, 0) return Arg1; } -DEFINE_PRIMITIVE ("IMAGE-SET-ROW!", Prim_image_set_row, 3, 3, 0) -{ long Length, i,j; - Pointer Pnrows, Pncols, Prest, Parray; - long nrows, ncols, row_to_set; - REAL *Array, *Row_Array; - - Primitive_3_Args(); - Arg_1_Type(TC_LIST); /* image = (nrows ncols array) */ - Pnrows = Vector_Ref(Arg1, CONS_CAR); - Prest = Vector_Ref(Arg1, CONS_CDR); - Pncols = Vector_Ref(Prest, CONS_CAR); - Prest = Vector_Ref(Prest, CONS_CDR); - Parray = Vector_Ref(Prest, CONS_CAR); - if (Vector_Ref(Prest, CONS_CDR) != NIL) Primitive_Error(ERR_ARG_1_WRONG_TYPE); - if (Type_Code(Parray) != TC_ARRAY) Primitive_Error(ERR_ARG_1_WRONG_TYPE); - - Range_Check(nrows, Pnrows, 0, 1024, ERR_ARG_1_BAD_RANGE); - Range_Check(ncols, Pncols, 0, 1024, ERR_ARG_1_BAD_RANGE); - Arg_2_Type(TC_FIXNUM); - Range_Check(row_to_set, Arg2, 0, (nrows-1), ERR_ARG_2_BAD_RANGE); - Arg_3_Type(TC_ARRAY); - Row_Array = Scheme_Array_To_C_Array(Arg3); - if (Array_Length(Arg3)>ncols) Primitive_Error(ERR_ARG_3_BAD_RANGE); - - Array = Scheme_Array_To_C_Array(Parray); - C_Image_Set_Row(Array, row_to_set, Row_Array, nrows, ncols); - return Arg1; + + +DEFINE_PRIMITIVE ("SUBIMAGE-COPY!", + Prim_subimage_copy, 12,12, 0) +{ long r1,c1, r2,c2, at1r,at1c,at2r,at2c, mr,mc; + long nn; + REAL *x,*y; + void subimage_copy(); + PRIMITIVE_HEADER (12); + CHECK_ARG (1, FIXNUM_P); /* rows 1 */ + CHECK_ARG (2, FIXNUM_P); /* cols 1 */ + CHECK_ARG (3, ARRAY_P); /* image array 1 = source array */ + CHECK_ARG (4, FIXNUM_P); /* rows 2 */ + CHECK_ARG (5, FIXNUM_P); /* cols 2 */ + CHECK_ARG (6, ARRAY_P); /* image array 2 = destination array */ + + CHECK_ARG (7, FIXNUM_P); /* at1 row */ + CHECK_ARG (8, FIXNUM_P); /* at1 col */ + CHECK_ARG (9, FIXNUM_P); /* at2 row */ + CHECK_ARG (10, FIXNUM_P); /* at2 col */ + CHECK_ARG (11, FIXNUM_P); /* m row */ + CHECK_ARG (12, FIXNUM_P); /* m col */ + + x = Scheme_Array_To_C_Array(ARG_REF(3)); + y = Scheme_Array_To_C_Array(ARG_REF(6)); + r1 = arg_nonnegative_integer(1); + c1 = arg_nonnegative_integer(2); + r2 = arg_nonnegative_integer(4); + c2 = arg_nonnegative_integer(5); + + nn = r1*c1; + if (nn != Array_Length(ARG_REF(3))) error_bad_range_arg(3); + nn = r2*c2; + if (nn != Array_Length(ARG_REF(6))) error_bad_range_arg(6); + + at1r = arg_nonnegative_integer(7); + at1c = arg_nonnegative_integer(8); + at2r = arg_nonnegative_integer(9); + at2c = arg_nonnegative_integer(10); + mr = arg_nonnegative_integer(11); + mc = arg_nonnegative_integer(12); + + if (((at1r+mr)>r1) || ((at1c+mc)>c1)) error_bad_range_arg(7); + if (((at2r+mr)>r2) || ((at2c+mc)>c2)) error_bad_range_arg(9); + + subimage_copy(x,y, r1,c1,r2,c2, at1r,at1c,at2r,at2c, mr,mc); + + PRIMITIVE_RETURN (NIL); } -DEFINE_PRIMITIVE ("IMAGE-SET-COLUMN!", Prim_image_set_column, 3, 3, 0) -{ long Length, i,j; - Pointer Pnrows, Pncols, Prest, Parray; - long nrows, ncols, col_to_set; - REAL *Array, *Col_Array; - - Primitive_3_Args(); - Arg_1_Type(TC_LIST); /* image = (nrows ncols array) */ - Pnrows = Vector_Ref(Arg1, CONS_CAR); - Prest = Vector_Ref(Arg1, CONS_CDR); - Pncols = Vector_Ref(Prest, CONS_CAR); - Prest = Vector_Ref(Prest, CONS_CDR); - Parray = Vector_Ref(Prest, CONS_CAR); - if (Vector_Ref(Prest, CONS_CDR) != NIL) Primitive_Error(ERR_ARG_1_WRONG_TYPE); - if (Type_Code(Parray) != TC_ARRAY) Primitive_Error(ERR_ARG_1_WRONG_TYPE); - - Range_Check(nrows, Pnrows, 0, 1024, ERR_ARG_1_BAD_RANGE); - Range_Check(ncols, Pncols, 0, 1024, ERR_ARG_1_BAD_RANGE); +void subimage_copy(x,y, r1,c1,r2,c2, at1r,at1c,at2r,at2c, mr,mc) + REAL *x,*y; long r1,c1,r2,c2, at1r,at1c,at2r,at2c, mr,mc; +{ long i,j; + REAL *xrow,*yrow; - Arg_2_Type(TC_FIXNUM); - Range_Check(col_to_set, Arg2, 0, (nrows-1), ERR_ARG_2_BAD_RANGE); - Arg_3_Type(TC_ARRAY); - Col_Array = Scheme_Array_To_C_Array(Arg3); - if (Array_Length(Arg3)>ncols) Primitive_Error(ERR_ARG_3_BAD_RANGE); + xrow = x + at1r*c1 + at1c; + yrow = y + at2r*c2 + at2c; /* A(i,j)--->Array[i*ncols+j] */ - Array = Scheme_Array_To_C_Array(Parray); - C_Image_Set_Col(Array, col_to_set, Col_Array, nrows, ncols); - return Arg1; + for (i=0; i Array[i*ncols + j] With no knowledge outside boundary, assume laplace(edge-point)=0.0 (no wrap-around, no artificial bndry) */ -C_image_laplacian(array, new_array, nrows, ncols) - REAL *array, *new_array; +void image_laplacian(x,y, nrows,ncols) + REAL *x, *y; long nrows, ncols; { long i,j, nrows1, ncols1; nrows1=nrows-1; ncols1=ncols-1; - if ((nrows<2)||(ncols<2)) return(1); /* no need todo anything for 1-point image */ + if ((nrows<2)||(ncols<2)) return; /* no need todo anything for 1-point image */ /* */ - i=0;j=0; new_array[i*ncols+j] = 0.0; /* NE corner */ - i=0;j=ncols1; new_array[i*ncols+j] = 0.0; /* NW corner */ - i=nrows1;j=0; new_array[i*ncols+j] = 0.0; /* SE corner */ - i=nrows1;j=ncols1; new_array[i*ncols+j] = 0.0; /* SW corner */ - i=0; for (j=1;j B(j,i) . @@ -986,7 +908,7 @@ Image_Transpose(Array, New_Array, nrows, ncols) New_Array[j*nrows + i] = Array[i*ncols + j]; /* (columns transposed-image) = nrows */ }} } - + /* IMAGE_ROTATE_90CLW A(i,j) <-> A(j, (nrows-1)-i) . @@ -1002,7 +924,7 @@ Image_Rotate_90clw(Array, Rotated_Array, nrows, ncols) Rotated_Array[(j*nrows) + ((nrows-1)-i)] = Array[i*ncols+j]; /* (columns rotated_image) =nrows */ }} } - + /* ROTATION 90degrees COUNTER-CLOCK-WISE: A(i,j) <-> A((nrows-1)-j, i) . (minus 1 because we start from 0). @@ -1021,7 +943,7 @@ Image_Rotate_90cclw(Array, Rotated_Array, nrows, ncols) Rotated_Array[to_index] = Array[from_index]; }} } - + /* IMAGE_MIRROR: A(i,j) <-> A(i, (ncols-1)-j) [ The -1 is there because we count from 0] . @@ -1043,8 +965,6 @@ C_Mirror_Image(Array, nrows, ncols) REAL *Array; long nrows, ncols; }} } - - /* IMAGE_ROTATE_90CLW_MIRROR: A(i,j) <-> A(j, i) this should be identical to image_transpose (see above). @@ -1063,384 +983,456 @@ C_Rotate_90clw_Mirror_Image(Array, Rotated_Array, nrows, ncols) Rotated_Array[to] = Array[from]; }} } - - - -/* END */ +/* More Image Manipulation ----------------------- + */ +DEFINE_PRIMITIVE ("SQUARE-IMAGE-TIME-REVERSE!", + Prim_square_image_time_reverse, 2,2, 0) +{ long i, rows; + REAL *a; + void square_image_time_reverse(); + PRIMITIVE_HEADER (2); + CHECK_ARG (1, ARRAY_P); + CHECK_ARG (2, FIXNUM_P); + a = Scheme_Array_To_C_Array(ARG_REF(1)); + rows = arg_nonnegative_integer(2); + if ((rows*rows) != Array_Length(ARG_REF(1))) error_bad_range_arg(1); + square_image_time_reverse(a,rows); + + PRIMITIVE_RETURN (NIL); +} +/* Square Image Time reverse + is combination of one-dimensional time-reverse + row-wise and column-wise. + It can be done slightly more efficiently than below. + */ +void square_image_time_reverse(x,rows) + REAL *x; + long rows; +{ long i,cols; + REAL *xrow, *yrow; + void C_Array_Time_Reverse(); + cols = rows; /* square image */ + + xrow = x; + for (i=0; i0.0)) return(.08 + .46 * (1 - t_bar)); - else return (0); -} - -REAL hanning(t, length) REAL t, length; -{ REAL twopi = 6.28318530717958; - REAL pi = twopi/2.; - REAL t_bar = cos(twopi * (t / length)); - if ((t0.0)) - return(.5 * (1 - t_bar)); - else return (0); -} - -REAL unit_square_wave(t) REAL t; -{ REAL twopi = 6.28318530717958; - REAL fmod(), fabs(); - REAL pi = twopi/2.; - REAL t_bar = fabs(fmod(t, twopi)); - if (t_bar < pi) return(1); - else return(0); + +/* From now on (below), assume that cs-images (rows=cols) have always EVEN LENGTH + which is true when they come from FFTs + */ + + + +/* In the following 3 time-reverse the bottom half rows + is done to match the frequencies of complex-images + coming from cft2d. + Also transpose is needed to match frequencies identically + + #| + ;; Scrabling of frequencies in cs-images + + ;; start from real image 4x4 + + ;; rft2d is a cs-image + (3.5 .375 -2.75 1.875 -.25 0. 0. -.25 -.25 -.125 0. .125 .25 .25 0. 0.) + + ;; cft2d transposed + ;; real + 3.5 .375 -2.75 .375 + -.25 0. 0. -.25 ; same as cs-image + -.25 -.125 0. -.125 + -.25 -.25 0. 0. ; row3 = copy 1 + time-reverse + ;; imag + 0. 1.875 0. -1.875 + .25 .25 0. 0. ; same as cs-image + 0. .125 0. -.125 + -.25 0. 0. -.25 ; row 3 = copy 1 + negate + time-reverse + |# + + */ + +void cs_image_magnitude(x,rows) + REAL *x; + long rows; +{ long i,j, cols, n,n2, nj; /* result = real ordinary image */ + REAL *xrow, *yrow; + cols = rows; /* input cs-image is square */ + n = rows; + n2 = n/2; + + xrow = x; + cs_array_magnitude(xrow, n); /* row 0 is cs-array */ + xrow = x + n2*cols; + cs_array_magnitude(xrow, n); /* row n2 is cs-array */ + + xrow = x + cols; /* real part */ + yrow = x + (rows-1)*cols; /* imag part */ + for (i=1; i