# HG changeset patch
# User jwe
# Date 775174745 0
# Node ID 4e826edfbc5662926d44d8db427710bc66f770e2
# Parent  faf108b99d21df4a19cbe55d1d80c9970cdf0667
[project @ 1994-07-25 22:18:28 by jwe]
Initial revision

diff --git a/scripts/general/prepad.m b/scripts/general/prepad.m
new file mode 100644
--- /dev/null
+++ b/scripts/general/prepad.m
@@ -0,0 +1,43 @@
+function y = prepad(x,l,c)
+#prepad(x,l)
+#Prepends zeros to the vector x until it is of length l.
+#prepad(x,l,c) prepends the constant c instead of zero.
+#
+#If length(x) > l, elements from the beginning of x are removed
+#until a vector of length l is obtained.
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+
+  if(nargin == 2)
+    c = 0;
+  elseif(nargin<2 || nargin>3)
+    error("usage: prepad(x,l) or prepad(x,l,c)");
+  endif
+
+  if(is_matrix(x))
+    error("first argument must be a vector");
+  elseif(!is_scalar(l))
+    error("second argument must be a scaler");
+  endif
+
+  if(l<0)
+    error("second argument must be non-negative");
+  endif
+
+  lx = length(x);
+
+  if(lx >= l)
+    y = x(lx-l+1:lx);
+  else
+    if(rows(x)>1)
+      y = [ c*ones(l-lx,1); x ];
+    else
+      y = [ c*ones(1,l-lx) x ];
+    endif
+  endif
+
+endfunction
diff --git a/scripts/image/colormap.m b/scripts/image/colormap.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/colormap.m
@@ -0,0 +1,41 @@
+function cmap = colormap(map)
+#Set the current colormap.
+#
+#colormap(map) sets the current colormap to map.  map should be an n row
+#by 3 column matrix. The columns contain red, green, and blue intensities
+#respectively.  All entries should be between 0 and 1 inclusive. The new
+#colormap is returned.
+#
+#colormap("default") restores the default colormap (a gray scale colormap
+#with 64 entries). The default colormap is returned.
+#
+#colormap with no arguments returns the current colormap.
+
+#Author:
+# Tony Richardson
+# amr@mpl.ucsd.edu
+# July 1994
+
+  global CURRENT_COLOR_MAP
+
+  cmap_name = "CURRENT_COLOR_MAP";
+
+  if(nargin == 1)
+    if(isstr(map))
+      if(strcmp(map,"default"))
+        CURRENT_COLOR_MAP = gray;
+      else
+        error("invalid argument");
+      endif
+    else
+      # Set the new color map
+      CURRENT_COLOR_MAP = map;
+    endif
+  elseif(exist(cmap_name) == 0)
+    # If global color map doesn't exist, create the default map.
+    CURRENT_COLOR_MAP = gray;
+  endif
+
+  # Return current color map.
+  cmap = CURRENT_COLOR_MAP;
+endfunction
diff --git a/scripts/image/gray.m b/scripts/image/gray.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/gray.m
@@ -0,0 +1,12 @@
+function map = gray(number)
+
+  if(nargin == 0)
+    number = 64;
+  endif
+
+  gr = [0:(number-1)]';
+  
+
+  map = [ gr gr gr ]/(number-1);
+
+endfunction
diff --git a/scripts/image/gray2ind.m b/scripts/image/gray2ind.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/gray2ind.m
@@ -0,0 +1,11 @@
+function [X, map] = gray2ind(I,n)
+
+  if(nargin == 1)
+    n = 64;
+  endif
+
+  map = gray(n);
+
+  X = round(I*(n-1)) + 1;
+
+endfunction
diff --git a/scripts/image/image.m b/scripts/image/image.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/image.m
@@ -0,0 +1,53 @@
+function image(x, zoom)
+#Display an octave image matrix.
+#
+#image(x) displays a matrix as a color image. The elements of x are indices
+#into the current colormap and should have values between 1 and the length
+#of the colormap.
+#
+#image(x,zoom) changes the zoom factor.  The default value is 4.
+#
+#SEE ALSO: imshow, imagesc, colormap.
+
+#Author:
+# Tony Richardson
+# amr@mpl.ucsd.edu
+# July 1994
+#
+#Modified:
+# Tony Richardson
+# amr@mpl.ucsd.edu
+# July 1994
+# (Modifications based on suggestions from John Eaton.)
+
+  global IMAGEDIR
+ 
+  if (nargin == 0)
+    # Load Bobbie Jo Richardson (Born 3/16/94)
+    x = loadimage([IMAGEDIR,"/default.img"]);
+    zoom = 2;
+  elseif(nargin == 1)
+    zoom = 4;
+  elseif(nargin > 2)
+    error("usage: image (matrix, [zoom])");
+  endif
+
+  # Generate random file name
+  rnd_str = num2str(fix(rand*10000));
+  ppm_name = ["image.", rnd_str, ".ppm" ];
+
+  saveimage(ppm_name,x,"ppm");
+
+  # Start the viewer
+  # Try xv, then xloadimage.
+
+  xv = sprintf("xv -expand %f %s",zoom,ppm_name);
+  xloadimage = sprintf("xloadimage -zoom %f %s",zoom*100, ppm_name);
+  rm = sprintf("rm -f %s",ppm_name);
+
+  command = sprintf("( %s || %s && %s ) > /dev/null 2>&1 &", ...
+                    xv, xloadimage, rm);
+
+  shell_cmd(command);
+
+endfunction
diff --git a/scripts/image/imagesc.m b/scripts/image/imagesc.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/imagesc.m
@@ -0,0 +1,47 @@
+function x = imagesc(x, zoom)
+#Scale and display a matrix as an image.
+#
+#imagesc(x) displays a scaled version of the matrix x.  The matrix is
+#scaled so that its entries are indices into the current colormap.
+#The scaled matrix is returned.
+#
+#imagesc(x,zoom) sets the magnification, the default value is 4.
+#
+#SEE ALSO: image, imshow
+
+#Author:
+# Tony Richardson
+# amr@mpl.ucsd.edu
+# July 1994
+#
+#Modified:
+# Tony Richardson
+# amr@mpl.ucsd.edu
+# July 1994
+# (Modifications based on suggestions from John Eaton.)
+
+  if (nargin < 1 || nargin > 2)
+    error("usage: image (matrix, [zoom])");
+  endif
+
+  if (nargin == 1)
+    zoom = 4;
+  endif
+
+  [ high, wide ] = size(x);
+
+  maxval = max(max(x));
+  minval = min(min(x));
+
+  # Rescale matrix so that all values are in the range 0 to
+  # length(colormap) inclusive
+  if (maxval == minval)
+    x = ones(high, wide);
+  else
+    # Rescale values to between 1 and length(colormap) inclusive.
+    x = fix((x - minval)/(maxval - minval) * (length(colormap)-1)) + 1;
+  endif
+
+  image(x,zoom);
+
+endfunction
diff --git a/scripts/image/imshow.m b/scripts/image/imshow.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/imshow.m
@@ -0,0 +1,33 @@
+function imshow(a1,a2,a3)
+#Display images.
+#
+#imshow(X) displays an indexed image using the current colormap.
+#
+#imshow(X,map) displays an indexed image using the specified colormap.
+#
+#imshow(I,n) displays a gray scale intensity image.
+#
+#imshow(R,G,B) displays an RGB image.
+#
+#SEE ALSO: image, imagesc, colormap, gray2ind, rgb2ind.
+
+#Author:
+# Tony Richardson
+# amr@mpl.ucsd.edu
+# July 1994
+
+  if (nargin == 0)
+    error("usage: imshow requires at least one argument.");
+  elseif(nargin == 2)
+    if(length(a2)==1)
+      [a1 a2] = gray2ind(a1,a2);
+    endif
+    colormap(a2);
+  elseif(nargin == 3)
+    [a1 a2] = rgb2ind(a1,a2,a3);
+    colormap(a2);
+  endif
+
+  image(a1);
+
+endfunction
diff --git a/scripts/image/ind2gray.m b/scripts/image/ind2gray.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/ind2gray.m
@@ -0,0 +1,35 @@
+function Y = ind2gray(X,map)
+#Convert an octave indexed image to a gray scale intensity image.
+#
+#Y = ind2gray(X) converts an indexed image to a gray scale intensity
+#image.  The current colormap is used to determine the intensities.
+#The intensity values lie between 0 and 1 inclusive.
+#
+#Y = ind2gray(X,map) uses the specified colormap instead of the current
+#one in the conversion process.
+#
+#SEE ALSO: gray2ind, rgb2ntsc, image, colormap
+
+  if (nargin == 1)
+    map = colormap;
+  endif
+
+  # Convert colormap to intensity values.
+  yiq = rgb2ntsc(map);
+  y = yiq(:,1);
+
+  # We need Fortran indexing capability, but be sure to save the user's
+  # preference.
+  pref = do_fortran_indexing;
+  do_fortran_indexing = "true";
+
+  # Replace indices in the input matrix with indexed values in the output
+  # matrix.
+  [rows, cols] = size(X);
+  Y = y(X(:));
+  Y = reshape(Y,rows,cols);
+
+  # Restore the user's preference.
+  do_fortran_indexing = pref;
+
+endfunction
diff --git a/scripts/image/ind2rgb.m b/scripts/image/ind2rgb.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/ind2rgb.m
@@ -0,0 +1,29 @@
+function [R G B] = ind2rgb(X,map)
+#Convert an indexed image to red, green, and blue color components.
+#
+#[R G B] = ind2rgb(X) uses the current colormap for the conversion.
+#
+#[R G B] = ind2rgb(X,map) uses the specified colormap.
+#
+#SEE ALSO: rgb2ind, image, imshow, ind2gray, gray2ind.
+
+  if(nargin == 1)
+    map = colormap;
+  endif
+
+  [hi wi] = size(X);
+
+  pref = do_fortran_indexing;
+  do_fortran_indexing = "true";
+
+  R = map(X(:),1);
+  G = map(X(:),2);
+  B = map(X(:),3);
+
+  R = reshape(R,hi,wi);
+  G = reshape(G,hi,wi);
+  B = reshape(B,hi,wi);
+
+  do_fortran_indexing = pref;
+
+endfunction
diff --git a/scripts/image/loadimage.m b/scripts/image/loadimage.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/loadimage.m
@@ -0,0 +1,12 @@
+function [X, map] = loadimage(filename)
+#Load an image file.
+#
+#[X, map] = loadimage(img_file) loads an image and it's associated color
+#map from file img_file.  The image must be in stored in octave's image
+#format.
+#
+#SEE ALSO: saveimage, load, save
+
+  eval(['load ', filename]);
+
+endfunction
diff --git a/scripts/image/ntsc2rgb.m b/scripts/image/ntsc2rgb.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/ntsc2rgb.m
@@ -0,0 +1,7 @@
+function rgb = ntsc2rgb(yiq)
+
+  trans = [ 1. 1. 1.; 0.95617 -0.27269 -1.10374; 0.62143 -0.64681 1.70062 ];
+
+  rgb = yiq * trans;
+
+endfunction
diff --git a/scripts/image/ocean.m b/scripts/image/ocean.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/ocean.m
@@ -0,0 +1,18 @@
+function map = ocean(number)
+
+  if(nargin == 0)
+    number = 64;
+  endif
+
+  cutin = fix(number/3);
+
+  
+  dr = (number-1)/(cutin);
+  r = prepad([0:dr:(number-1)],number)';
+  dg = (number-1)/(2*cutin);
+  g = prepad([0:dg:(number-1)],number)';
+  b = [0:(number-1)]';
+
+  map = [ r g b ]/(number-1);
+
+endfunction
diff --git a/scripts/image/octtopnm.c b/scripts/image/octtopnm.c
new file mode 100644
--- /dev/null
+++ b/scripts/image/octtopnm.c
@@ -0,0 +1,302 @@
+# /*
+cc -s -o octtopnm octtopnm.c
+exit
+*/
+
+#include <stdio.h>
+#include <string.h>
+#include <malloc.h>
+
+/* usage: octtopnm [-a] octfile */
+
+static void usage(message)
+char *message;
+{
+  if(message != NULL) {
+    fprintf(stderr,"octtopnm: %s\n",message);
+  }
+  fprintf(stderr,"usage: octtopnm [-a] octavefile\n");
+  exit(1);
+}
+
+static void fatal(message)
+char *message;
+{
+  if(message != NULL) {
+    fprintf(stderr,"octtopnm: %s\n",message);
+  }
+  exit(1);
+}
+
+int main(argc, argv)
+int argc;
+char **argv;
+{
+  int rawbits = 1, row, col, index;
+  int cmap_rows, cmap_cols, img_rows, img_cols;
+  int gray, pbm, pgm, ppm;
+  unsigned char **rgb, byte;
+  unsigned short **img;
+  char *oct_file_name;
+  FILE *oct_file;
+  char cmap_name[4], cmap_type[7], img_name[2], img_type[7];
+  double mat_val;
+  int option;
+  extern char *optarg;
+  extern int optind;
+
+  if(argc == 1) {
+    usage(NULL);
+  }
+
+  while((option = getopt(argc,argv,"ha")) != EOF) {
+    switch(option) {
+    case 'h':
+      /* help */
+      usage(NULL);
+      break;
+    case 'a':
+      rawbits = 0;
+      break;
+    case '?':
+    default:
+      usage("unrecognized option");
+    }
+  }
+
+  if(optind+1 != argc) {
+    usage("input file name missing");
+  }
+
+  oct_file_name = argv[optind];
+  if((oct_file = fopen(oct_file_name,"r")) == NULL) {
+    fatal("unable to open input file");
+  }
+
+  if(fscanf(oct_file,"# name: %s\n",cmap_name) != 1 || 
+     strcmp(cmap_name,"map") != 0) {
+    fatal("not a valid octave image file");
+  }
+
+  if(fscanf(oct_file,"# type: %s\n",cmap_type) != 1 || 
+     strcmp(cmap_type,"matrix") != 0) {
+    fatal("not a valid octave image file");
+  }
+
+  if(fscanf(oct_file,"# rows: %d\n",&cmap_rows) != 1) {
+    fatal("error reading octave image file");
+  }
+
+  if(fscanf(oct_file,"# columns: %d\n",&cmap_cols) != 1) {
+    fatal("error reading octave image file");
+  }
+
+  if(cmap_cols != 3) {
+    fatal("invalid color map in octave image file");
+  }
+
+  if((rgb = (unsigned char **)
+      malloc(cmap_rows*sizeof(unsigned char *))) == NULL) {
+    fatal("out of memory");
+  }
+
+  if((rgb[0] = (unsigned char *)
+      malloc(cmap_rows*cmap_cols*sizeof(unsigned char))) == NULL) {
+    fatal("out of memory");
+  }
+
+  for(row=1; row<cmap_rows; row++) {
+    rgb[row] = rgb[row-1]+3;
+  }
+
+  gray = 1;
+  for(row=0; row<cmap_rows; row++) {
+    for(col=0; col<cmap_cols; col++) {
+      if(fscanf(oct_file,"%lf",&mat_val) != 1) {
+        fatal("error reading color map entries");
+      }
+      if(mat_val < 0) mat_val = 0.;
+      if(mat_val > 1) mat_val = 1.;
+      rgb[row][col] = mat_val*255;
+    }
+    if(gray) {
+      if(rgb[row][0] != rgb[row][1] || rgb[row][0] != rgb[row][2]) {
+        /* It's a color image. */
+        gray = 0;
+      }
+    }
+  }
+
+  if(fscanf(oct_file,"\n# name: %s\n",img_name) != 1 || 
+     strcmp(img_name,"X") != 0) {
+    fatal("not a valid octave image file");
+  }
+
+  if(fscanf(oct_file,"# type: %s\n",img_type) != 1 || 
+     strcmp(img_type,"matrix") != 0) {
+    fatal("not a valid octave image file");
+  }
+
+  if(fscanf(oct_file,"# rows: %d\n",&img_rows) != 1) {
+    fatal("error reading octave image file");
+  }
+
+  if(fscanf(oct_file,"# columns: %d\n",&img_cols) != 1) {
+    fatal("error reading octave image file");
+  }
+
+  if((img = (unsigned short **)
+      malloc(img_rows*sizeof(unsigned short *))) == NULL) {
+    fatal("out of memory");
+  }
+
+  if((img[0] = (unsigned short *)
+      malloc(img_rows*img_cols*sizeof(unsigned short))) == NULL) {
+    fatal("out of memory");
+  }
+
+  for(row=1; row<img_rows; row++) {
+    img[row] = img[row-1]+img_cols;
+  }
+
+  for(row=0; row<img_rows; row++) {
+    for(col=0; col<img_cols; col++) {
+      if(fscanf(oct_file,"%lf",&mat_val) != 1) {
+        fatal("error reading color map entries");
+      }
+      if(mat_val < 1) mat_val = 1.;
+      if(mat_val > cmap_rows) mat_val = cmap_rows;
+      img[row][col] = mat_val;
+    }
+  }
+
+  pbm = pgm = ppm = 0;
+
+  if(cmap_rows == 2 && gray && 
+     ((rgb[0][0] == 0 && rgb[1][0] == 255) ||
+      (rgb[0][0] == 255 && rgb[1][0] == 0))) {
+    /* Create a bitmap only if there are two colormap entries and they are
+       black and white. */
+    pbm = 1;
+  }
+  else if(gray) {
+    /* If not a bitmap, create a gray scale image if the entries within
+       each row of the color map are equal. */
+    pgm = 1;
+  }
+  else {
+    /* Otherwise create a full color image. */
+    ppm = 1;
+  }
+
+  if(rawbits) {
+    if(pbm) {
+      printf("P4\n");
+      printf("%d %d\n",img_cols,img_rows);
+      index = 0;
+      for(row=0; row<img_rows; row++) {
+        for(col=0; col<img_cols; col++) {
+          if(index == 7) {
+            byte = (byte << 1) + !rgb[img[row][col]-1][0];
+            fwrite(&byte,sizeof(unsigned char),1,stdout);
+            byte = 0;
+            index = 0;
+          }
+          else {
+            byte = (byte << 1) + !rgb[img[row][col]-1][0];
+            index++;
+          }
+        }
+      }
+      if(index != 0) {
+        printf("\n");
+      }
+    }
+    else if(pgm) {
+      printf("P5\n");
+      printf("%d %d\n",img_cols,img_rows);
+      printf("255\n");
+      for(row=0; row<img_rows; row++) {
+        for(col=0; col<img_cols; col++) {
+          fwrite(rgb[img[row][col]-1],sizeof(unsigned char),1,stdout);
+        }
+      }
+    }
+    else {
+      printf("P6\n");
+      printf("%d %d\n",img_cols,img_rows);
+      printf("255\n");
+      for(row=0; row<img_rows; row++) {
+        for(col=0; col<img_cols; col++) {
+          fwrite(rgb[img[row][col]-1],sizeof(unsigned char),3,stdout);
+        }
+      }
+    }
+  }
+  else {
+    if(pbm) {
+      printf("P1\n");
+      printf("%d %d\n",img_cols,img_rows);
+      index = 0;
+      for(row=0; row<img_rows; row++) {
+        for(col=0; col<img_cols; col++) {
+          if(index == 30) {
+            printf("%d\n",!rgb[img[row][col]-1][0]);
+            index = 0;
+          }
+          else {
+            printf("%d ",!rgb[img[row][col]-1][0]);
+            index++;
+          }
+        }
+      }
+      if(index != 0) {
+        printf("\n");
+      }
+    }
+    else if(pgm) {
+      printf("P2\n");
+      printf("%d %d\n",img_cols,img_rows);
+      printf("255\n");
+      index = 0;
+      for(row=0; row<img_rows; row++) {
+        for(col=0; col<img_cols; col++) {
+          if(index == 12) {
+            printf("%d\n",rgb[img[row][col]-1][0]);
+            index = 0;
+          }
+          else {
+            printf("%d ",rgb[img[row][col]-1][0]);
+            index++;
+          }
+        }
+      }
+      if(index != 0) {
+        printf("\n");
+      }
+    }
+    else {
+      printf("P3\n");
+      printf("%d %d\n",img_cols,img_rows);
+      printf("255\n");
+      index = 0;
+      for(row=0; row<img_rows; row++) {
+        for(col=0; col<img_cols; col++) {
+          if(index == 4) {
+            printf("%d %d %d\n",rgb[img[row][col]-1][0],
+                   rgb[img[row][col]-1][1],rgb[img[row][col]-1][2]);
+            index = 0;
+          }
+          else {
+            printf("%d %d %d ",rgb[img[row][col]-1][0],
+                   rgb[img[row][col]-1][1],rgb[img[row][col]-1][2]);
+            index++;
+          }
+        }
+      }
+      if(index != 0) {
+        printf("\n");
+      }
+    }
+  }
+}
diff --git a/scripts/image/rgb2ind.m b/scripts/image/rgb2ind.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/rgb2ind.m
@@ -0,0 +1,31 @@
+function [X map] = rgb2ind(R,G,B)
+#Convert and RGB image to an octave indexed image.
+#
+#[X map] = rgb2ind(R,G,B)
+#
+#SEE ALSO: ind2rgb, rgb2ntsc.
+#
+#Bugs: The color map may have duplicate entries.
+
+  if(nargin != 3)
+    error("usage: [X, map] = rgb2ind(R,G,B)");
+  endif
+
+  [hi wi] = size(R);
+
+  X = zeros(hi,wi);
+
+  map = zeros(hi*wi,3);
+
+  pref = do_fortran_indexing;
+  do_fortran_indexing = "true";
+
+  map(:,1) = R(:);
+  map(:,2) = G(:);
+  map(:,3) = B(:);
+
+  X(:) = 1:(hi*wi);
+
+  do_fortran_indexing = pref;
+
+endfunction
diff --git a/scripts/image/rgb2ntsc.m b/scripts/image/rgb2ntsc.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/rgb2ntsc.m
@@ -0,0 +1,6 @@
+function yiq = rgb2ntsc(rgb)
+
+  trans = [ 0.299 0.596 0.211; 0.587 -0.274 -0.523; 0.114 -0.322 0.312 ];
+  yiq = rgb * trans;
+
+endfunction
diff --git a/scripts/image/saveimage.m b/scripts/image/saveimage.m
new file mode 100644
--- /dev/null
+++ b/scripts/image/saveimage.m
@@ -0,0 +1,75 @@
+function saveimage(filename,X,img_form,map)
+#Save a matrix to disk in image format.
+#
+#saveimage(filename,x) saves matrix x to file filename in octave's image
+#format.  The current colormap is saved in the file also.
+#
+#saveimage(filename,x,"img") saves the image in the default format and
+#is the same as saveimage(filename,x).
+#
+#saveimage(filename,x,"ppm") saves the image in ppm format instead of
+#the default octave image format.
+#
+#saveimage(filename,x,"ps") saves the image in PostScript format instead
+#of the default octave image format. (Note: images saved in PostScript format
+#can not be read back into octave with loadimage.)
+#
+#saveimage(filename,x,format,map) saves the image along with the specified
+#colormap in the specified format.
+#
+#Note: If the colormap contains only two entries and these entries are black
+#and white, the bitmap ppm and PostScript formats are used.  If the image is
+#a gray scale image (the entries within each row of the colormap are equal)
+#the gray scale ppm and PostScript image formats are used, otherwise the full
+#color formats are used.
+#
+#SEE ALSO: loadimage, save, load, colormap
+
+  global IMAGEDIR
+
+  if(nargin < 2)
+    error("usage: saveimage(filename,matrix,[format, [colormap]])");
+  elseif(nargin == 2)
+    if(!isstr(filename))
+      error("File name must be a string.");
+    endif
+    map = colormap;
+    img_form = "img";
+  elseif(nargin == 3)
+    if(!isstr(img_form))
+      error("Image format specification must be a string");
+    endif
+    map = colormap;
+  endif
+
+  if (strcmp(img_form,"img") == 1)
+    oct_file = filename;
+  elseif (strcmp(img_form,"ppm") == 1)
+    # We need a temporary octave image file name.
+    oct_file = sprintf("image.%s.img",num2str(fix(rand*10000)));
+    ppm_file = filename;
+  elseif (strcmp(img_form,"ps") == 1)
+    # We need a temporary octave image file name.
+    oct_file = sprintf("image.%s.img",num2str(fix(rand*10000)));
+    ps_file = filename;
+  endif
+
+  # Save image in octave image file format
+  eval(['save ', oct_file, ' map X']);
+
+  # Convert to another format if requested.
+  if (strcmp(img_form,"ppm") == 1)
+    octtopnm = sprintf([IMAGEDIR,"/octtopnm %s > %s"],oct_file,filename);
+    rm = sprintf("rm -f %s",oct_file);
+    shell_cmd(octtopnm);
+    shell_cmd(rm);
+  elseif (strcmp(img_form,"ps") == 1)
+    octtopnm = sprintf([IMAGEDIR,"/octtopnm %s"],oct_file);
+    ppmtops = sprintf("pnmtops > %s 2> /dev/null", filename);
+    octtops = [ octtopnm, " | ", ppmtops ];
+    rm = sprintf("rm -f %s",oct_file);
+    shell_cmd(octtops);
+    shell_cmd(rm);
+  endif
+
+endfunction
diff --git a/scripts/polynomial/polyinteg.m b/scripts/polynomial/polyinteg.m
new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/polyinteg.m
@@ -0,0 +1,38 @@
+function p = polyinteg(p)
+#polyinteg(c)
+#Returns the coefficients of the integral the polynomial whose coefficients
+#are represented by the vector c.
+#
+#The constant of integration is zero.
+#
+#SEE ALSO: poly, polyderiv, polyreduce, roots, conv, deconv, residue,
+#          filter, polyval, polyvalm
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+  if(nargin != 1)
+    error("usage: polyinteg(vector)");
+  endif
+
+  if(is_matrix(p))
+    error("argument must be a vector");
+  endif
+
+  lp = length(p);
+
+  if(lp == 0)
+    p = [];
+    return;
+  end
+
+  if(rows(p) > 1)
+    # Convert to column vector
+    p = p.';
+  endif
+
+  p = [ p 0 ] ./ [lp:-1:1 1];
+
+endfunction
diff --git a/scripts/polynomial/residue.m b/scripts/polynomial/residue.m
new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/residue.m
@@ -0,0 +1,268 @@
+function [r, p, k, e] = residue(b,a,toler)
+#[r p k e] = residue(b,a)
+#If b and a are vectors of polynomial coefficients, then residue
+#calculates the partial fraction expansion corresponding to the
+#ratio of the two polynomials. The vector r contains the residue
+#terms, p contains the pole values, k contains the coefficients of
+#a direct polynomial term (if it exists) and e is a vector containing
+#the powers of the denominators in the partial fraction terms.
+#Assuming b and a represent polynomials P(s) and Q(s) we have:
+#
+# P(s)    M       r(m)         N
+# ---- =  #  -------------  +  # k(n)*s^(N-n)
+# Q(s)   m=1 (s-p(m))^e(m)    n=1
+#
+#(# represents summation) where M is the number of poles (the length of
+#the r, p, and e vectors) and N is the length of the k vector.
+#
+#[r p k e] = residue(b,a,tol)
+#This form of the function call may be used to set a tolerance value.
+#The default value is 0.001. The tolerance value is used to determine
+#whether poles with small imaginary components are declared real. It is
+#also used to determine if two poles are distinct. If the ratio of the
+#imaginary part of a pole to the real part is less than tol, the imaginary
+#part is discarded. If two poles are farther apart than tol they are
+#distinct.
+#
+#Example:
+# b = [1 1 1];
+# a = [1  -5   8  -4];
+#
+# [r p k e] = residue(b,a)
+#
+# returns
+#
+# r = [-2 7 3]; p = [2 2 1]; k = []; e = [1 2 1];
+#
+# which implies the following partial fraction expansion
+#
+#       s^2 + s + 1         -2       7        3
+#   ------------------- = ----- + ------- + -----
+#   s^3 - 5s^2 + 8s - 4   (s-2)   (s-2)^2   (s-1)
+#
+#SEE ALSO: poly, roots, conv, deconv, polyval, polyderiv, polyinteg
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+# Here's the method used to find the residues.
+# The partial fraction expansion can be written as:
+#
+#          
+#   P(s)    D   M(k)      A(k,m)
+#   ---- =  #    #    -------------
+#   Q(s)   k=1  m=1   (s - pr(k))^m
+#
+# (# is used to represent a summation) where D is the number of
+# distinct roots, pr(k) is the kth distinct root, M(k) is the
+# multiplicity of the root, and A(k,m) is the residue cooresponding
+# to the kth distinct root with multiplicity m.  For example,
+#
+#        s^2            A(1,1)   A(2,1)    A(2,2)
+# ------------------- = ------ + ------ + -------
+# s^3 + 4s^2 + 5s + 2    (s+2)    (s+1)   (s+1)^2
+#
+# In this case there are two distinct roots (D=2 and pr = [-2 -1]),
+# the first root has multiplicity one and the second multiplicity
+# two (M = [1 2]) The residues are actually stored in vector format as
+# r = [ A(1,1) A(2,1) A(2,2) ].
+#
+# We then multiply both sides by Q(s).  Continuing the example:
+#
+# s^2 = r(1)*(s+1)^2 + r(2)*(s+1)*(s+2) + r(3)*(s+2)
+#
+# or
+#
+# s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2)
+#
+# The coefficients of the polynomials on the right are stored
+# in a row vector called rhs, while the coefficients of the
+# polynomial on the left is stored in a row vector called lhs.
+# If the multiplicity of any root is greater than one we'll
+# also need derivatives of this equation of order up to the
+# maximum multiplicity minus one.  The derivative coefficients
+# are stored in successive rows of lhs and rhs.
+#
+# For our example lhs and rhs would be:
+#
+#       | 1 0 0 |
+# lhs = |       |
+#       | 0 2 0 |
+#
+#       | 1 2 1 1 3 2 0 1 2 |
+# rhs = |                   |
+#       | 0 2 2 0 2 3 0 0 1 |
+#
+# We then form a vector B and a matrix A obtained by evaluating the
+# polynomials in lhs and rhs at the pole values. If a pole has a
+# multiplicity greater than one we also evaluate the derivative
+# polynomials (successive rows) at the pole value.
+#
+# For our example we would have
+#
+# | 4|   | 1 0 0 |   | r(1) |
+# | 1| = | 0 0 1 | * | r(2) |
+# |-2|   | 0 1 1 |   | r(3) |
+#
+# We then solve for the residues using matrix division.
+
+  if(nargin < 2 || nargin > 3)
+    error("usage: residue(b,a[,toler])");
+  endif
+
+  if (nargin == 2)
+    # Set the default tolerance level
+    toler = .001;
+  endif
+
+  # Make sure both polynomials are in reduced form.
+  a = polyreduce(a);
+  b = polyreduce(b);
+
+  b = b/a(1);
+  a = a/a(1);
+
+  la = length(a);
+  lb = length(b);
+
+  # Handle special cases here.
+  if(la == 0 || lb == 0)
+    k = r = p = e = [];
+    return;
+  elseif (la == 1)
+    k = b/a;
+    r = p = e = [];
+    return;
+  endif
+
+  # Find the poles.
+  p = roots(a);
+  lp = length(p);
+
+  # Determine if the poles are (effectively) real.
+  index = find(abs(imag(p) ./ real(p)) < toler);
+  if (length(index) != 0)
+    p(index) = real(p(index));
+  endif
+
+  # Find the direct term if there is one.
+  if(lb>=la)
+    # Also returns the reduced numerator.
+    [k, b] = deconv(b,a);
+    lb = length(b);
+  else
+    k = [];
+  endif
+
+  if(lp == 1)
+    r = polyval(b,p);
+    e = 1;
+    return;
+  endif
+
+
+  # We need to determine the number and multiplicity of the roots.
+  # D is the number of distinct roots.
+  # M is a vector of length D containing the multiplicity of each root.
+  # pr is a vector of length D containing only the distinct roots.
+  # e is a vector of length lp which indicates the power in the partial
+  # fraction expansion of each term in p.
+
+  # Set initial values.  We'll remove elements from pr as we find
+  # multiplicities.  We'll shorten M afterwards.
+  e = ones(lp,1);
+  M = zeros(lp,1);
+  pr = p;
+  D = 1; M(1) = 1;
+
+  old_p_index = 1; new_p_index = 2;
+  M_index = 1; pr_index = 2;
+  while(new_p_index<=lp)
+    if(abs(p(new_p_index) - p(old_p_index)) < toler)
+      # We've found a multiple pole.
+      M(M_index) = M(M_index) + 1;
+      e(new_p_index) = e(new_p_index-1) + 1;
+      # Remove the pole from pr.
+      pr(pr_index) = [];
+    else
+      # It's a different pole.
+      D++; M_index++; M(M_index) = 1;
+      old_p_index = new_p_index; pr_index++;
+    endif
+    new_p_index++;
+  endwhile
+
+  # Shorten M to it's proper length
+  M = M(1:D);
+
+  # Now set up the polynomial matrices.
+
+  MM = max(M);
+  # Left hand side polynomial
+  lhs = zeros(MM,lb);
+  rhs = zeros(MM,lp*lp);
+  lhs(1,:) = b;
+  rhi = 1;
+  dpi = 1;
+  mpi = 1;
+  while(dpi<=D)
+    for ind = 1:M(dpi)
+      if(mpi>1 && (mpi+ind)<=lp)
+        cp = [p(1:mpi-1); p(mpi+ind:lp)];
+      elseif (mpi==1)
+        cp = p(mpi+ind:lp);
+      else
+        cp = p(1:mpi-1);
+      endif
+      rhs(1,rhi:rhi+lp-1) = prepad(poly(cp),lp);
+      rhi = rhi + lp;
+    endfor
+    mpi = mpi + M(dpi);
+    dpi++;
+  endwhile
+  if(MM > 1)
+    for index = 2:MM
+      lhs(index,:) = prepad(polyderiv(lhs(index-1,:)),lb);
+      ind = 1;
+      for rhi = 1:lp
+        cp = rhs(index-1,ind:ind+lp-1);
+        rhs(index,ind:ind+lp-1) = prepad(polyderiv(cp),lp);
+        ind = ind + lp;
+      endfor
+    endfor
+  endif
+
+  # Now lhs contains the numerator polynomial and as many derivatives as are
+  # required.  rhs is a matrix of polynomials, the first row contains the
+  # corresponding polynomial for each residue and successive rows are
+  # derivatives.
+
+  # Now we need to evaluate the first row of lhs and rhs at each distinct
+  # pole value.  If there are multiple poles we will also need to evaluate
+  # the derivatives at the pole value also.
+
+  B = zeros(lp,1);
+  A = zeros(lp,lp);
+
+  dpi = 1;
+  row = 1;
+  while(dpi<=D)
+    for mi = 1:M(dpi)
+      B(row) = polyval(lhs(mi,:),pr(dpi));
+      ci = 1;
+      for col = 1:lp
+        cp = rhs(mi,ci:ci+lp-1);
+        A(row,col) = polyval(cp,pr(dpi));
+        ci = ci + lp;
+      endfor
+      row++;
+    endfor
+    dpi++;
+  endwhile
+
+  # Solve for the residues.
+  r = A\B;
+
+endfunction
diff --git a/scripts/polynomial/roots-amr.m b/scripts/polynomial/roots-amr.m
new file mode 100644
--- /dev/null
+++ b/scripts/polynomial/roots-amr.m
@@ -0,0 +1,55 @@
+function r = roots(c)
+#Find the roots of a polynomial.
+#
+#In octave, a polynomial is represented by it's coefficients (arranged
+#in descending order). For example, a vector c of length n+1 corresponds
+#to the following nth order polynomial
+#
+#  p(x) = c(1) x^n + ... + c(n) x + c(n+1).
+#
+#roots(c) will return a vector r of length n such that
+#
+#  p(x) = c(1) [ (x - r(1)) * (x - r(2)) * ... * (x - r(n)) ]
+#
+#roots and poly are inverse functions to within a scaling factor.
+#
+#SEE ALSO: poly, roots, conv, deconv, residue, filter, polyval, polyvalm,
+#          polyderiv, polyinteg
+
+# Author:
+#  Tony Richardson
+#  amr@mpl.ucsd.edu
+#  June 1994
+
+  if(nargin != 1)
+    error("usage: roots(vector)");
+  endif
+
+  if(is_matrix(c))
+    error("argument must be a vector.");
+  endif
+
+  n = length(c);
+
+  if(is_scalar(c) || n == 0)
+    r = [];
+    return;
+  endif
+
+  # Ensure that c is a row vector.
+  if(rows(c) > 1)
+    c = reshape(c,1,n);
+  endif
+
+  # We could replace this with a call to compan, but it's faster to
+  # just reproduce the code here.
+  A = diag(ones(n-2,1),-1);
+  A(1,:) = -c(2:n)/c(1);
+
+  r = eig(A);
+ 
+  # Sort roots in order by decreasing magnitude.
+  [mr i] = sort(abs(r));
+  r = r(i(length(i):-1:1));
+
+endfunction
diff --git a/scripts/set/complement.m b/scripts/set/complement.m
new file mode 100644
--- /dev/null
+++ b/scripts/set/complement.m
@@ -0,0 +1,31 @@
+function y = complement(a,b)
+
+# usage: complement(a,b)
+#
+# Returns the elements of set b that are not in set a.
+#
+# See - create_set, union, intersection
+
+  if(nargin != 2)
+    error("usage: complement(a,b)");
+  endif
+
+  if(isempty(a))
+    y = create_set(b);
+  elseif(isempty(b))
+    y = [];
+  else
+    a = create_set(a);
+    b = create_set(b);
+    yindex = 1;
+    y = zeros(1,length(b));
+    for index = 1:length(b)
+      if(all(a != b(index)))
+        y(yindex++) = b(index);
+      endif
+    endfor
+    y = y(1:(yindex-1));
+  endif
+
+endfunction
+
diff --git a/scripts/set/create_set.m b/scripts/set/create_set.m
new file mode 100644
--- /dev/null
+++ b/scripts/set/create_set.m
@@ -0,0 +1,36 @@
+function y = create_set(x)
+
+# usage: create_set(x)
+#
+# Returns the unique elements of x, sorted in ascending order.
+#
+# See - union, intersection, complement
+
+  if ( nargin != 1)
+    error("usage: create_set(x)");
+  endif
+
+  if(isempty(x))
+    y = [];
+  else
+    [nrx, ncx] = size(x);
+    nelx = nrx*ncx;
+    x = reshape(x,1,nelx);
+    y = zeros(1,nelx);
+
+    x = sort(x);
+    cur_val = y(1) = x(1);
+    yindex = xindex = 2;
+
+    while (xindex <= nelx)
+      if(cur_val != x(xindex))
+        cur_val = x(xindex);
+        y(yindex++) = cur_val;
+      endif
+      xindex++;
+    endwhile
+    y = y(1:(yindex-1));
+  endif
+  
+endfunction
+  
diff --git a/scripts/set/intersection.m b/scripts/set/intersection.m
new file mode 100644
--- /dev/null
+++ b/scripts/set/intersection.m
@@ -0,0 +1,41 @@
+function y = intersection(a,b)
+
+# usage: intersection(a,b)
+#
+# Returns the intersection of sets a and b.
+#
+# See - create_set, union, complement
+
+  if (nargin != 2)
+    error("usage: intersection(a,b)");
+  endif
+
+  if(isempty(a) || isempty(b))
+    y = [];
+    return;
+  endif
+
+  a = create_set(a);
+  b = create_set(b);
+
+  if(length(a) < length(b))
+    yindex = 1;
+    y = zeros(1,length(a));
+    for index = 1:length(a)
+      if(any(b == a(index)))
+        y(yindex++) = a(index);
+      endif
+    endfor
+  else
+    yindex = 1;
+    y = zeros(1,length(b));
+    for index = 1:length(b)
+      if(any(a == b(index)))
+        y(yindex++) = b(index);
+      endif
+    endfor
+  endif
+
+  y = y(1:(yindex-1));
+
+endfunction
diff --git a/scripts/set/union.m b/scripts/set/union.m
new file mode 100644
--- /dev/null
+++ b/scripts/set/union.m
@@ -0,0 +1,25 @@
+function y = union(a,b)
+
+# usage: union(a,b)
+#
+# Returns the union of sets a and b.
+#
+# See - create_set, intersection, complement
+
+  if (nargin != 2)
+    error("usage: union(a,b)");
+  endif
+
+  if(isempty(a))
+    y = create_set(b);
+  elseif(isempty(b))
+    y = create_set(a);
+  else
+    [nra, nca] = size(a);
+    a = reshape(a,1,nra*nca);
+    [nrb, ncb] = size(b);
+    b = reshape(b,1,nrb*ncb);
+    y = create_set([a, b]);
+  endif
+
+endfunction
diff --git a/scripts/signal/sinc.m b/scripts/signal/sinc.m
new file mode 100644
--- /dev/null
+++ b/scripts/signal/sinc.m
@@ -0,0 +1,31 @@
+function result = sinc(x)
+
+# usage: sinc(x)
+#
+#        Returns sin(pi*x)/(pi*x).
+#
+
+# We either need to set the do_fortran_indexing variable to "true"
+# or use reshape to convert the input matrix to a vector, so that
+# we can use find to determine the elements of x that equal zero.
+# I prefer reshaping.
+
+  [nr, nc] = size(x);
+
+  nels = nr*nc;
+
+  x = reshape(x,nels,1);
+
+  # Set result to all ones initially.
+  result = ones(nels,1);
+
+  # Find non-zero elements in the input matrix.
+  i = find(x);
+
+  if (!isempty(i))
+    result(i) = sin(pi*x(i))./(pi*x(i));
+  endif
+
+  result = reshape(result,nr,nc);
+
+endfunction