unit geomprnt;

interface

uses
  arrtypes, Math, SysUtils, Windows;


procedure i4mat_transpose_print(m: Integer; n: Integer; var a: TIntegerArray; title: String);

procedure i4mat_transpose_print_some(m: Integer; n: Integer; var a: TIntegerArray; ilo: Integer; jlo: Integer; ihi: Integer; jhi: Integer; title: string);

procedure r8mat_transpose_print(m: Integer; n: Integer; a: TRealArray; title: string);

procedure r8mat_transpose_print_some(m: Integer; n: Integer; a: TRealArray; ilo: Integer; jlo: Integer; ihi: Integer; jhi: Integer; title: string);

procedure r8vec_print(n: Integer; var a: TRealArray; title: String);

procedure timestamp();

function timestring(): string;

function triangulation_plot_eps(file_out_name: string; g_num: Integer; g_xy: TRealArray; tri_num: Integer; nod_tri: TIntegerArray): Boolean;

procedure triangulation_print(point_num: Integer; xc: TRealArray; tri_num: Integer; tri_vert: TIntegerArray; tri_nabe: TIntegerArray);



implementation

uses
  geompack;


procedure i4mat_transpose_print(m: Integer; n: Integer; var a: TIntegerArray; title: String);

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_TRANSPOSE_PRINT prints an I4MAT, transposed.
//
//  Modified:
//
//    31 January 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows in A.
//
//    Input, int N, the number of columns in A.
//
//    Input, int A[M*N], the M by N matrix.
//
//    Input, char *TITLE, a title to be printed.
//
begin
  i4mat_transpose_print_some(m, n, a, 1, 1, m, n, title);
end;
//****************************************************************************80

procedure i4mat_transpose_print_some(m: Integer; n: Integer; var a: TIntegerArray; ilo: Integer; jlo: Integer; ihi: Integer; jhi: Integer; title: string);

//****************************************************************************80
//
//  Purpose:
//
//    I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed.
//
//  Modified:
//
//    09 February 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows of the matrix.
//    M must be positive.
//
//    Input, int N, the number of columns of the matrix.
//    N must be positive.
//
//    Input, int A[M*N], the matrix.
//
//    Input, int ILO, JLO, IHI, JHI, designate the first row and
//    column, and the last row and column to be printed.
//
//    Input, char *TITLE, a title for the matrix.

var
  i: Integer;
  i2hi: Integer;
  i2lo: Integer;
  j: Integer;
  j2hi: Integer;
  j2lo: Integer;
const
  INCX: Integer = 10;
begin
  if (0 < Length(Trim(title))) then begin
    WriteLn('');
    WriteLn(title);
  end;
//
//  Print the columns of the matrix, in strips of INCX.
//
  i2lo:= ilo;
  while(i2lo <= ihi) do begin
    i2hi:= i2lo + INCX - 1;
    i2hi:= Min(i2hi, m);
    i2hi:= Min(i2hi, ihi);

    WriteLn('');
//
//  For each row I in the current range...
//
//  Write the header.
//
    Write('  Row:    ');
    for i:= i2lo to i2hi do begin
      Write(Format('%7d       ',[i]));
    end;
    WriteLn('');
    WriteLn('  Col');
    WriteLn('');
//
//  Determine the range of the rows in this strip.
//
    j2lo:= Max(jlo, 1);
    j2hi:= Min(jhi, n);

    for j:= j2lo to j2hi do begin
//
//  Print out(up to INCX) entries in column J, that lie in the current strip.
//
      Write(Format('%5d  ',[j]));
      for i:= i2lo to i2hi do begin
        Write(Format('%6d  ',[a[i-1+(j-1)*m]]));
      end;
      WriteLn('');
    end;
    inc(i2lo, INCX);
  end;

  WriteLn('');
end;
//****************************************************************************80

procedure r8mat_transpose_print(m: Integer; n: Integer; a: TRealArray; title: string);

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_TRANSPOSE_PRINT prints a real matrix, transposed.
//
//  Modified:
//
//    11 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Input, double A[M*N], an M by N matrix to be printed.
//
//    Input, char *TITLE, an optional title.
//
begin
  r8mat_transpose_print_some(m, n, a, 1, 1, m, n, title);
end;
//****************************************************************************80

procedure r8mat_transpose_print_some(m: Integer; n: Integer; a: TRealArray; ilo: Integer; jlo: Integer; ihi: Integer; jhi: Integer; title: string);

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_TRANSPOSE_PRINT_SOME prints some of a real matrix, transposed.
//
//  Modified:
//
//    11 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Input, double A[M*N], an M by N matrix to be printed.
//
//    Input, int ILO, JLO, the first row and column to print.
//
//    Input, int IHI, JHI, the last row and column to print.
//
//    Input, char *TITLE, an optional title.
//
var
  i: Integer;
  i2: Integer;
  i2hi: Integer;
  i2lo: Integer;
  inc: Integer;
  j: Integer;
  j2hi: Integer;
  j2lo: Integer;
const
  INCX: Integer = 5;
begin
  if (0 < Length(Trim(title))) then begin
    WriteLn('');
    WriteLn(title);
  end;


  i2lo:= Max(ilo, 1);
  while (i2lo <= Min(ihi, m)) do begin
    i2hi:= i2lo + INCX - 1;
    i2hi:= Min(i2hi, m);
    i2hi:= Min(i2hi, ihi);

    inc:= i2hi + 1 - i2lo;

    WriteLn('');
    Write('  Row: ');
    for i:= i2lo to i2hi do begin
      Write(Format('%7d       ',[i]));
    end;
    WriteLn('');
    WriteLn('  Col');
    WriteLn('');

    j2lo:= Max(jlo, 1);
    j2hi:= Min(jhi, n);

    for j:= j2lo to j2hi do begin
      Write(Format('%5d ',[j]));
      for i2:= 1 to inc do begin
        i:= i2lo - 1 + i2;
        Write(Format('%14.5f',[a[(i-1)+(j-1)*m]]));
      end;
      WriteLn('');
    end;

    i2lo:= i2lo + INCX;
  end;

  WriteLn('');
end;
//****************************************************************************80

procedure r8vec_print(n: Integer; var a: TRealArray; title: String);

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_PRINT prints an R8VEC.
//
//  Modified:
//
//    16 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of components of the vector.
//
//    Input, double A[N], the vector to be printed.
//
//    Input, char *TITLE, a title to be printed first.
//    TITLE may be blank.
//
var
  i: Integer;
begin
  if (0 < Length(Trim(title))) then begin
    WriteLn('');
    WriteLn(title);
  end;

  WriteLn('');
  for i:= 0 to n-1 do begin
    WriteLn(Format('%6d   %14.5f',[i + 1, a[i]]));
  end;
end;
//****************************************************************************80

procedure timestamp();

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Modified:
//
//    21 August 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
var
  f: TFormatSettings;
  loc: Cardinal;
begin
  loc:= (SORT_DEFAULT shl 16) or (SUBLANG_ENGLISH_US shl 10) or LANG_ENGLISH;
  GetLocaleFormatSettings(loc, f);
  f.ShortDateFormat:= 'dd MMMM yyyy';
  f.LongTimeFormat:= 'hh:mm:ss AMPM';

  WriteLn(DateTimeToStr(Now(), f));
end;
//****************************************************************************80

function timestring(): string;

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTRING returns the current YMDHMS date as a string.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Modified:
//
//    13 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *TIMESTRING, a string containing the current YMDHMS date.
//
var
  f: TFormatSettings;
  loc: Cardinal;
begin
  loc:= (SORT_DEFAULT shl 16) or (SUBLANG_ENGLISH_US shl 10) or LANG_ENGLISH;
  GetLocaleFormatSettings(loc, f);
  f.ShortDateFormat:= 'dd MMMM yyyy';
  f.LongTimeFormat:= 'hh:mm:ss AMPM';

  Result:= DateTimeToStr(Now(), f);
end;
//****************************************************************************80

function triangulation_plot_eps(file_out_name: string; g_num: Integer; g_xy: TRealArray; tri_num: Integer; nod_tri: TIntegerArray): Boolean;

//****************************************************************************80
//
//  Purpose:
//
//    TRIANGULATION_PLOT_EPS plots a triangulation of a pointset.
//
//  Discussion:
//
//    The triangulation is most usually a Delaunay triangulation,
//    but this is not necessary.
//
//    The data can be generated by calling DTRIS2, but this is not
//    necessary.
//
//  Modified:
//
//    08 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *FILE_OUT_NAME, the name of the output file.
//
//    Input, int G_NUM, the number of points.
//
//    Input, double G_XY[G_NUM,2], the coordinates of the points.
//
//    Input, int TRI_NUM, the number of triangles.
//
//    Input, int NOD_TRI[3,TRI_NUM], lists, for each triangle,
//    the indices of the points that form the vertices of the triangle.
//
//    Output, bool TRIANGULATION_PLOT_EPS, is TRUE for success.
//
var
  date_time: string;
  e: Integer;
  file_out: TextFile;
  g: Integer;
  j: Integer;
  k: Integer;
  t: Integer;
  x_max: Real;
  x_min: Real;
  x_ps: Integer;
  x_ps_max: Integer;
  x_ps_max_clip: Integer;
  x_ps_min: Integer;
  x_ps_min_clip: Integer;
  y_max: Real;
  y_min: Real;
  y_ps: Integer;
  y_ps_max: Integer;
  y_ps_max_clip: Integer;
  y_ps_min: Integer;
  y_ps_min_clip: Integer;
begin
  x_ps_max:= 576;
  x_ps_max_clip:= 594;
  x_ps_min:= 36;
  x_ps_min_clip:= 18;
  y_ps_max:= 666;
  y_ps_max_clip:= 684;
  y_ps_min:= 126;
  y_ps_min_clip:= 108;

  date_time:= timestring();

  x_max:= g_xy[0+0*2];
  x_min:= g_xy[0+0*2];
  y_max:= g_xy[1+0*2];
  y_min:= g_xy[1+0*2];

  for g:= 0 to g_num-1 do begin
    x_max:= Max(x_max, g_xy[0+g*2]);
    x_min:= Min(x_min, g_xy[0+g*2]);
    y_max:= Max(y_max, g_xy[1+g*2]);
    y_min:= Min(y_min, g_xy[1+g*2]);
  end;
//
//  Plot the Delaunay triangulation.
//
//
//  Open the output file.
//
  AssignFile(file_out, file_out_name);
  Rewrite(file_out);

  WriteLn(file_out, '%!PS-Adobe-3.0 EPSF-3.0');
  WriteLn(file_out, '%%Creator: triangulation_plot_eps.cc');
  WriteLn(file_out, '%%Title: ' + file_out_name);
  WriteLn(file_out, '%%CreationDate: '+ date_time);
  WriteLn(file_out, '%%Pages: 1');
  WriteLn(file_out, '%%Bounding Box: ' +IntToStr(x_ps_min) + '  ' + IntToStr(y_ps_min) + '  ' + IntToStr(x_ps_max) + '  ' + IntToStr(y_ps_max));
  WriteLn(file_out, '%%Document-Fonts: Times-Roman');
  WriteLn(file_out, '%%LanguageLevel: 1');
  WriteLn(file_out, '%%EndComments');
  WriteLn(file_out, '%%BeginProlog');
  WriteLn(file_out, '/inch {72 mul} def');
  WriteLn(file_out, '%%EndProlog');
  WriteLn(file_out, '%%Page: 1 1');
  WriteLn(file_out, 'save');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%  Set the RGB line color to very light gray.');
  WriteLn(file_out, '%');
  WriteLn(file_out, '0.900  0.900  0.900 setrgbcolor');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%  Draw a gray border around the page.');
  WriteLn(file_out, '%');
  WriteLn(file_out, 'newpath');
  WriteLn(file_out, '  ' + IntToStr(x_ps_min) + '  ' + IntToStr(y_ps_min) + ' moveto');
  WriteLn(file_out, '  ' + IntToStr(x_ps_max) + '  ' + IntToStr(y_ps_min) + ' lineto');
  WriteLn(file_out, '  ' + IntToStr(x_ps_max) + '  ' + IntToStr(y_ps_max) + ' lineto');
  WriteLn(file_out, '  ' + IntToStr(x_ps_min) + '  ' + IntToStr(y_ps_max) + ' lineto');
  WriteLn(file_out, '  ' + IntToStr(x_ps_min) + '  ' + IntToStr(y_ps_min) + ' lineto');
  WriteLn(file_out, 'stroke');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%  Set the RGB line color to black.');
  WriteLn(file_out, '%');
  WriteLn(file_out, '0.000  0.000  0.000 setrgbcolor');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%  Set the font and its size.');
  WriteLn(file_out, '%');
  WriteLn(file_out, '/Times-Roman findfont');
  WriteLn(file_out, '0.50 inch scalefont');
  WriteLn(file_out, 'setfont');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%  Print a title.');
  WriteLn(file_out, '%');
  WriteLn(file_out, '210  702  moveto');
  WriteLn(file_out, '(Triangulation)  show');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%  Define a clipping polygon.');
  WriteLn(file_out, '%');
  WriteLn(file_out, 'newpath');
  WriteLn(file_out, '  ' + IntToStr(x_ps_min_clip) + '  ' + IntToStr(y_ps_min_clip) + ' moveto');
  WriteLn(file_out, '  ' + IntToStr(x_ps_max_clip) + '  ' + IntToStr(y_ps_min_clip) + ' lineto');
  WriteLn(file_out, '  ' + IntToStr(x_ps_max_clip) + '  ' + IntToStr(y_ps_max_clip) + ' lineto');
  WriteLn(file_out, '  ' + IntToStr(x_ps_min_clip) + '  ' + IntToStr(y_ps_max_clip) + ' lineto');
  WriteLn(file_out, '  ' + IntToStr(x_ps_min_clip) + '  ' + IntToStr(y_ps_min_clip) + ' lineto');
  WriteLn(file_out, 'clip newpath');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%  Set the RGB line color to green.');
  WriteLn(file_out, '%');
  WriteLn(file_out, '0.000  0.750  0.150 setrgbcolor');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%  Draw the nodes.');
  WriteLn(file_out, '%');

  for g:= 0 to g_num-1 do begin
    x_ps:= Trunc(
      ( ( x_max - g_xy[0+g*2] ) * ( x_ps_min )
      + ( g_xy[0+g*2] - x_min ) * ( x_ps_max ) )
      / ( x_max - x_min ) );

    y_ps:= Trunc(
      ( ( y_max - g_xy[1+g*2] ) * ( y_ps_min )
      + ( g_xy[1+g*2] - y_min ) * ( y_ps_max ) )
      / ( y_max - y_min ) );

    WriteLn(file_out, 'newpath ' + IntToStr(x_ps) + '  '
             + IntToStr(y_ps) + ' 5 0 360 arc closepath fill');
  end;

  WriteLn(file_out, '%');
  WriteLn(file_out, '%  Set the RGB line color to red.');
  WriteLn(file_out, '%');
  WriteLn(file_out, '0.900  0.200  0.100 setrgbcolor');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%  Draw the triangles.');
  WriteLn(file_out, '%');

  for t:= 1 to tri_num do begin
    WriteLn(file_out, 'newpath');

    for j:= 1 to 4 do begin
      e:= i4_wrap(j, 1, 3);

      k:= nod_tri[3*(t-1)+e-1];

      x_ps:= Trunc(
        ( ( x_max - g_xy[0+(k-1)*2] ) * ( x_ps_min )
        + ( g_xy[0+(k-1)*2] - x_min ) * ( x_ps_max ) )
        / ( x_max - x_min ) );

      y_ps:= Trunc(
        ( ( y_max - g_xy[1+(k-1)*2] ) * ( y_ps_min )
        + ( g_xy[1+(k-1)*2] - y_min ) * ( y_ps_max ) )
        / ( y_max - y_min ) );

      if (j = 1) then begin
        WriteLn(file_out, IntToStr(x_ps) + '  ' + IntToStr(y_ps) + ' moveto');
      end
      else begin
        WriteLn(file_out, IntToStr(x_ps) + '  ' + IntToStr(y_ps) + ' lineto');
      end;
    end;

    WriteLn(file_out, 'stroke');
  end;

  WriteLn(file_out, 'restore  showpage');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%  End of page.');
  WriteLn(file_out, '%');
  WriteLn(file_out, '%%Trailer');
  WriteLn(file_out, '%%EOF');

  CloseFile(file_out);

  Result:= true;
end;
//****************************************************************************80

procedure triangulation_print(point_num: Integer; xc: TRealArray; tri_num: Integer; tri_vert: TIntegerArray; tri_nabe: TIntegerArray);

//****************************************************************************80
//
//  Purpose:
//
//    TRIANGULATION_PRINT prints information defining a Delaunay triangulation.
//
//  Discussion:
//
//    Triangulations created by RTRIS include extra information encoded
//    in the negative values of TRI_NABE.
//
//    Because some of the nodes counted in POINT_NUM may not actually be
//    used in the triangulation, I needed to compute the true number
//    of vertices.  I added this calculation on 13 October 2001.
//
//    Ernest Fasse pointed out an error in the indexing of VERTEX_LIST,
//    which was corrected on 19 February 2004.
//
//  Modified:
//
//    19 February 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double XC[2*POINT_NUM], the point coordinates.
//
//    Input, int TRI_NUM, the number of triangles.
//
//    Input, int TRI_VERT[3*TRI_NUM], the nodes that make up the triangles.
//
//    Input, int TRI_NABE[3*TRI_NUM], the triangle neighbors on each side.
//    If there is no triangle neighbor on a particular side, the value of
//    TRI_NABE should be negative.  If the triangulation data was created by
//    DTRIS2, then there is more information encoded in the negative values.
//
var
  boundary_num: Integer;
  i: Integer;
  j: Integer;
  k: Integer;
  n1: Integer;
  n2: Integer;
  s: Integer;
  s1: Integer;
  s2: Integer;
  skip: Boolean;
  t: Integer;
  vertex_list: TIntegerArray;
  vertex_num: Integer;
const
  DIM_NUM: Integer = 2;
begin
  WriteLn('');
  WriteLn('TRIANGULATION_PRINT');
  WriteLn('  Information defining a triangulation.');
  WriteLn('');
  WriteLn(Format('  The number of points is %d', [point_num]));

  r8mat_transpose_print(DIM_NUM, point_num, xc, '  Point coordinates');

  WriteLn('');
  WriteLn(Format('  The number of triangles is %d',[tri_num]));
  WriteLn('');
  WriteLn('  Sets of three points are used as vertices of');
  WriteLn('  the triangles.  For each triangle, the points');
  WriteLn('  are listed in counterclockwise order.');

  i4mat_transpose_print(3, tri_num, tri_vert, '  Triangle nodes');

  WriteLn('');
  WriteLn('  On each side of a given triangle, there is either');
  WriteLn('  another triangle, or a piece of the convex hull.');
  WriteLn('  For each triangle, we list the indices of the three');
  WriteLn('  neighbors, or(if negative) the codes of the');
  WriteLn('  segments of the convex hull.');

  i4mat_transpose_print(3, tri_num, tri_nabe, '  Triangle neighbors');
//
//  Determine VERTEX_NUM, the number of vertices.  This is not
//  the same as the number of points!
//
  SetLength(vertex_list, 3*tri_num);

  k:= 0;
  for t:= 0 to tri_num-1 do begin
    for s:= 0 to 2 do begin
      vertex_list[k]:= tri_vert[s+t*3];
      k:= k + 1;
    end;
  end;

  i4vec_sort_heap_a(3*tri_num, vertex_list);

  i4vec_sorted_unique(3*tri_num, vertex_list, vertex_num);
//
//  Determine the number of boundary points.
//
  boundary_num:= 2 * vertex_num - tri_num - 2;

  WriteLn('');
  WriteLn(Format('  The number of boundary points is %d', [boundary_num]));
  WriteLn('');
  WriteLn('  The segments that make up the convex hull can be');
  WriteLn('  determined from the negative entries of the triangle');
  WriteLn('  neighbor list.');
  WriteLn('');
  WriteLn('  #Tri Side  N1  N2');
  WriteLn('');

  skip:= false;

  k:= 0;

  for i:= 0 to tri_num-1 do begin
    for j:= 0 to 2 do begin
      if (tri_nabe[j+i*3] < 0) then begin
        s:= -tri_nabe[j+i*3];
        t:= s div 3;

        if ((t < 1) or (tri_num < t)) then begin
          WriteLn('');
          WriteLn('  Sorry, this data does not use the DTRIS2');
          WriteLn('  convention for convex hull segments.');
          skip:= true;
          break;
        end;

        s1:=(s mod 3) + 1;
        s2:= i4_wrap(s1+1, 1, 3);
        k:= k + 1;
        n1:= tri_vert[s1-1+(t-1)*3];
        n2:= tri_vert[s2-1+(t-1)*3];
        WriteLn(Format('%4d  %4d  %4d  %4d  %4d', [k,t,s1,n1,n2]));
      end;
    end;

    if (skip) then begin
      break;
    end;

  end;
end;

end.