octave_packages/vrml-1.0.13/vrml_surf.m

   1 ## Copyright (C) 2002-2009 Etienne Grossmann <etienne@egdn.net>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## s = vrml_surf (x, y, z [, options] ) - code for a VRML surface
  17 ## s = vrml_surf (z [, options] )
  18 ##
  19 ## Returns vrml97 code for a Shape -> IndexedFaceSet node representing a
  20 ## surface passing through the given points.
  21 ##
  22 ## x : RxC or C  : X coordinates of the points on the surface
  23 ## y : RxC or R  : Y "                                      "
  24 ## z : RxC       : Z "                                      "
  25 ##
  26 ## s :   string  : The code
  27 ##
  28 ## If x and y are omitted, they are assumed to be linspace(-1,1,C or R).
  29 ## Points presenting one or more 'inf' or 'nan' coordinates are ignored.
  30 ##
  31 ## Options :
  32 ##
  33 ## "col" , col  : 3      : RGB Color,                default = [0.3,0.4,0.9]
  34 ##             or 3x(R*C): Color of vertices (vrml colorPerVertex is TRUE).
  35 ##             or 3x((R-1)*(C-1))
  36 ##                       : Color of facets
  37 ##             or 1      : Reflectivity (equivalent to [col,col,col] in RGB)
  38 ##             or R x C  : Reflectivity of vertices
  39 ##             or 1x(R*C)
  40 ##             or (R-1)x(C-1)
  41 ##             or (R-1)*(C-1)
  42 ##                       : Reflectivity of facets.
  43 ##
  44 ##        RGB and reflectivity values should be in the [0,1] interval.
  45 ##
  46 ## "checker", c : 1x2 : Color as a checker. If c(1) is positive, checker has
  47 ##            c(1) rows. If it is negative, each checker row is c(1) facets
  48 ##            high c(2) likewise determines width of checker columns.
  49 ## "checker", c : 1x1 : Same as [c,c].
  50 ##
  51 ## "zcol", zc   : 3xN : Specify a colormap. The color of each vertex is
  52 ##            interpolated according to its height (z).
  53 ##
  54 ## "zgray"      : Black-to-white colormap. Same as "zcol", [0 1;0 1;0 1].
  55 ##
  56 ## "zrb"        : Red-to-blue. Same as "zcol", [0 7 10;0 0 2;7 19 2]/10.
  57 ##
  58 ## "steps"      : Represent surface as a piecewise constant Z = f(X,Y) function
  59 ##
  60 ## "bars"       : Represent surface as a bar plot
  61 ##
  62 ## "tran", tran : 1x1    : Transparency,                        default = 0
  63 ##
  64 ## "creaseAngle", a
  65 ##              : 1      : vrml creaseAngle The browser may smoothe the fold
  66 ##                         between facets forming an angle less than a.
  67 ##                                                              default = 0
  68 ## "smooth"           : same as "creaseAngle",pi.
  69 ## "tex", texFile
  70 ##
  71 ## See also: vmesh(), vrml_faces(), test_moving_surf()
  72
  73 function s = vrml_surf (x, y, z,varargin)
  74
  75 if (nargin <= 1) || ischar(y),  # Cruft to allow not passing x and y
  76   zz = x ;
  77   [R,C] = size (zz);
  78   [xx,yy] = meshgrid (linspace (-1,1,C), linspace (-1,1,R));
  79   ## ones(R,1)*[1:C] ;
  80   ## yy = ## [1:R]'*ones(1,C) ;
  81   ##if     nargin >=3,
  82   ##  s = vrml_surf ( xx, yy, zz, y, z, varargin{:} );
  83   ##  return
  84   ##elseif nargin >=2,
  85   ##  s = vrml_surf ( xx, yy, zz, y, varargin{:} );
  86   ##  return
  87   ##end
  88   if nargin >= 3
  89     varargin = {y, z, varargin{:}};
  90   elseif nargin >= 2
  91     varargin = {y, varargin{:}};
  92   end
  93
  94   x = xx ; y = yy ; z = zz ;
  95 end
  96
  97 defaultCol = [0.3; 0.4; 0.9];
  98
  99                                 # Read options
 100                                 # Default values
 101
 102 upper = 1;                      # Do "upper" triangulation of square grid
 103 tran = 0 ;                      # Transparency
 104 col = defaultCol ;              # Color
 105 checker = 0;                    # Checkered coloring
 106 colorPerVertex = 1;             # Color vertices or faces
 107 zrb = zgray = zcol = 0;         # Color by elevation
 108 emit = 0;                       # emissiveColor or diffuse only
 109 smooth = creaseAngle = nan ;
 110 steps = 0;
 111 bars = 0;
 112 tex = 0;
 113 bwid = -1;
 114 DEFcoord = DEFcol = "";         # Give a name to VRML objects
 115
 116 if numel (varargin)
 117
 118   op1 = " tran col creaseAngle emit colorPerVertex checker DEFcoord DEFcol zcol bwid tex ";
 119   op0 = " smooth zgray zrb steps bars " ;
 120
 121   default = tars (tran, col, creaseAngle, emit, colorPerVertex, steps, bars, \
 122                   bwid, DEFcoord, DEFcol, zcol, smooth, checker, zgray, zrb, tex);
 123
 124   s = read_options (varargin,"op0",op0,"op1",op1,"default",default);
 125
 126   tran=            s.tran;
 127   col=             s.col;
 128   creaseAngle=     s.creaseAngle;
 129   emit=            s.emit;
 130   colorPerVertex=  s.colorPerVertex;
 131   DEFcoord=        s.DEFcoord;
 132   DEFcol=          s.DEFcol;
 133   zcol=            s.zcol;
 134   smooth=          s.smooth;
 135   steps=           s.steps;
 136   bars=            s.bars;
 137   bwid=            s.bwid;
 138   tex=             s.tex;
 139   if bwid >= 0
 140     bars = 1;
 141   end
 142   checker=         s.checker;
 143   zgray=           s.zgray;
 144   zrb=             s.zrb;
 145 end
 146 ## keyboard
 147 if ! isnan (smooth), creaseAngle = pi ; end
 148
 149 [R,C] = size(z);
 150 if any (size (x) == 1), x = ones(R,1)*x(:)' ; end
 151 if any (size (y) == 1), y = y(:)*ones(1,C)  ; end
 152
 153 if bars
 154
 155   if  bwid < 0
 156     bwid = 2/3;
 157   end
 158   brad = bwid/2;
 159
 160   R4 = 4*R;
 161   C4 = 4*C;
 162   x2 = y2 = z2 = zeros (R4,C4);
 163
 164   x2(:,1) = x2(:,2) = kron ((1+brad)*x(:,1) - brad*x(:,2), [1;1;1;1]);
 165   x2(:,C4-1) = x2(:,C4) = kron ((1+brad)*x(:,C) - brad*x(:,C-1), [1;1;1;1]);
 166   x2(:,5:4:C4) =  x2(:,6:4:C4) = kron ((1-brad)*x(:,2:C) + brad*x(:,1:C-1), [1;1;1;1]);
 167   x2(:,3:4:C4-4) =  x2(:,4:4:C4-4) = kron ((1-brad)*x(:,1:C-1)+brad*x(:,2:C), [1;1;1;1]);
 168
 169   y2(1,:) = y2(2,:) = kron ((1+brad)*y(1,:) - brad*y(2,:), [1 1 1 1]);
 170   y2(R4-1,:) = y2(R4,:) = kron ((1+brad)*y(R,:) - brad*y(R-1,:), [1 1 1 1]);
 171   y2(5:4:R4,:) =  y2(6:4:R4,:) = kron ((1-brad)*y(2:R,:) + brad*y(1:R-1,:), [1 1 1 1]);
 172   y2(3:4:R4-4,:) =  y2(4:4:R4-4,:) = kron ((1-brad)*y(1:R-1,:) + brad*y(2:R,:), [1 1 1 1]);
 173
 174   z2([2:4:R4;3:4:R4],[2:4:C4;3:4:C4]) = kron(z,ones(2));
 175
 176   x = x2;
 177   y = y2;
 178   z = z2;
 179   R = R4;
 180   C = C4;
 181
 182   if numel (size (col)) == 2 && all (size (col) == size (defaultCol)) && all (col == defaultCol)
 183     col = [col, 0.8*col, 0.9*col];
 184   end
 185   if numel (col) == 3
 186     col = col(:);
 187     topCol = col;
 188     botCol = col;
 189     sideCol = [col,col,col,col];
 190   elseif numel (col) == 6
 191     col = col(:);
 192     topCol = col(1:3);
 193     botCol = col(1:3);
 194     sideCol = col(4:6)*ones(1,4);
 195   elseif numel (col) == 9
 196     col = col(:);
 197     topCol = col(1:3);
 198     botCol = col(7:9);
 199     sideCol = col(4:6)*ones(1,4);
 200   end
 201   col = ones(3, R-1, C-1);
 202   for i=1:3
 203     col(i,2:4:R-1,2:4:C-1) = topCol(i);
 204     col(i,4:4:R-1,:) = botCol(i);
 205     col(i,:,4:4:C-1) = botCol(i);
 206     col(i,3:4:R-1,2:4:C-1) = sideCol(i,1);
 207     col(i,2:4:R-1,1:4:C-1) = sideCol(i,2);
 208     col(i,1:4:R-1,2:4:C-1) = sideCol(i,3);
 209     col(i,2:4:R-1,3:4:C-1) = sideCol(i,2);
 210   end
 211   iOnFloor = find (! z(1:R-1,1:C-1));
 212   if ! isempty (iOnFloor)
 213     ## keyboard
 214     col(3*(iOnFloor-1)+1) = botCol(1);
 215     col(3*(iOnFloor-1)+2) = botCol(2);
 216     col(3*(iOnFloor-1)+3) = botCol(3);
 217   end
 218
 219 elseif steps                    # Constant by parts
 220
 221                                 # Intermediate coordinates  (R+1) x (C+1)
 222   x2 = (x([1,1:R],[1,1:C]) + x([1,1:R],[1:C,C])) / 2;
 223   y2 = (y([1,1:R],[1:C,C]) + y([1:R,R],[1:C,C])) / 2;
 224
 225                                 # Extend extremities so all patches have same size
 226   x2(1,:) = 2*x2(1,:) - x2(2,:);
 227   x2(:,1) = 2*x2(:,1) - x2(:,2);
 228   x2(R+1,:) = 2*x2(R+1,:) - x2(R,:);
 229   x2(:,C+1) = 2*x2(:,C+1) - x2(:,C);
 230
 231   y2(1,:) = 2*y2(1,:) - y2(2,:);
 232   y2(:,1) = 2*y2(:,1) - y2(:,2);
 233   y2(R+1,:) = 2*y2(R+1,:) - y2(R,:);
 234   y2(:,C+1) = 2*y2(:,C+1) - y2(:,C);
 235
 236                                 # Duplicate intermediate values 2R x 2C
 237   ii = [1,([1;1]*(2:R))(:)',R+1];
 238   jj = [1,([1;1]*(2:C))(:)',C+1];
 239
 240   x2 = x2(ii,jj);
 241   y2 = y2(ii,jj);
 242
 243   z2 = z([1;1]*(1:R),[1;1]*(1:C));;
 244   x = x2;
 245   y = y2;
 246   z = z2;
 247
 248   if checker
 249     col = checker_color (checker, col, 2*R,2*C);
 250   end
 251   if numel (col) == R*C
 252     col = [1;1;1]*col(:)';
 253   end
 254   if numel (col) == 3*R*C
 255     col = reshape (col, 3,R,C);
 256     col2 = zeros (3,2*R-1,2*C-1);
 257     col2(1,:,:) = defaultCol(1);
 258     col2(2,:,:) = defaultCol(2);
 259     col2(3,:,:) = defaultCol(3);
 260     col2(:,1:2:end,1:2:end) = col;
 261     col = reshape (col2,3,(2*R-1)*(2*C-1));
 262   end
 263
 264   R *= 2;
 265   C *= 2;
 266 end
 267
 268 pts = [x(:)';y(:)';z(:)'];
 269
 270 keepp = all (!isnan(pts) & finite(pts)) ;
 271 keepip = find (keepp);
 272 if tex
 273   [texX,texY] = meshgrid (linspace (0,1,C), linspace (0,1,R));
 274   texXY = [texX(:)'; texY(:)'];
 275 end
 276
 277 trgs = zeros(3,2*(R-1)*(C-1)) ;
 278
 279 ## Points are numbered as
 280 ##
 281 ## 1  R+1 .. (C-1)*R+1
 282 ## 2  R+2         :
 283 ## :   :          :
 284 ## R  2*R ..     C*R
 285 ##
 286
 287
 288 ## (x,y), (x,y+1), (x+1,y)  i.e. n, n+1, n+R
 289
 290 if !upper                       # Do regular triangulation
 291   ## Triangles are numbered as :
 292   ##
 293   ## X = (R-1)*(C-1)
 294   ## +-----+-----+-----+-----+-----+
 295   ## |    /|    /|    /|    /|-R+1/|
 296   ## | 1 / | R / |   / |   / |R*C/ |
 297   ## |  /  |  /  |  /  |  /  |  /  |
 298   ## | /X+1| /X+R| /   | /   | /   |
 299   ## |/    |/    |/    |/    |/    |
 300   ## +-----+-----+-----+-----+-----+
 301   ## |    /|    /|    /|    /|    /|
 302   ## | 2 / |R+2/ |   / |   / |   / |
 303   ## |  /  |  /  |  /  |  /  |  /  |
 304   ## | /   | /   | /   | /   | /   |
 305   ## |/    |/    |/    |/    |/    |
 306   ## +-----+-----+-----+-----+-----+
 307   ##    :           :           :
 308   ##    :           :           :
 309   ## +-----+-----+-----+-----+-----+
 310   ## |    /|    /|    /|    /|    /|
 311   ## |R-1/ |2*R/ |   / |   / |C*R/ |
 312   ## |  /  |-1/  |  /  |  /  |  /  |
 313   ## | /X+R| /   | /   | /   | /C*R|
 314   ## |/    |/    |/    |/    |/ X+ |
 315   ## +-----+-----+-----+-----+-----+
 316
 317   tmp = 1:(R-1)*(C-1);
 318
 319   trgs(1,tmp) = ([1:R-1]'*ones(1,C-1) + R*ones(R-1,1)*[0:C-2])(:)';
 320   trgs(2,tmp) = ([ 2:R ]'*ones(1,C-1) + R*ones(R-1,1)*[0:C-2])(:)';
 321   trgs(3,tmp) = ([1:R-1]'*ones(1,C-1) + R*ones(R-1,1)*[1:C-1])(:)';
 322
 323   tmp += (R-1)*(C-1);
 324
 325   trgs(1,tmp) = ([1:R-1]'*ones(1,C-1) + R*ones(R-1,1)*[1:C-1])(:)';
 326   trgs(2,tmp) = ([ 2:R ]'*ones(1,C-1) + R*ones(R-1,1)*[0:C-2])(:)';
 327   trgs(3,tmp) = ([ 2:R ]'*ones(1,C-1) + R*ones(R-1,1)*[1:C-1])(:)';
 328
 329 else                            # Do "upper" triangulation
 330
 331   ##  Each triangle is      +-----+     +-----+
 332   ##  the highest of either |    /| or  |\    |
 333   ##                        |   / |     | \   |
 334   ##                        |  /  |     |  \  |
 335   ##                        | /   |     |   \ |
 336   ##                        |/    |     |    \|
 337   ##                        +-----+     +-----+
 338
 339
 340   tmp = 1:(R-1)*(C-1);
 341   tmp2 = reshape(1:R*C,R,C);
 342   foo1 = z(1:R-1,1:C-1) + z(2:R,2:C);
 343   foo2 = z(2:R,1:C-1) + z(1:R-1,2:C);
 344   tmp3 = (!isnan(foo1) & (isnan (foo2) | foo1 > foo2))(:)';
 345
 346   trgs(1,tmp) = tmp2(1:R-1,1:C-1)(:)';
 347   trgs(2,tmp) = tmp2(2:R,1:C-1)(:)';
 348   trgs(3,tmp) = trgs(1,tmp) + R + tmp3 ;
 349   tmp += (R-1)*(C-1);
 350   trgs(1,tmp) = tmp2(1:R-1,2:C)(:)';
 351   trgs(2,tmp) = tmp2(2:R,2:C)(:)';
 352   trgs(3,tmp) = trgs(1,tmp) - R + 1 - tmp3 ;
 353 end                             # EOF "upper" triangulation
 354
 355 #trgs = trgs(:,find(rem(1:R-1,2)'*rem(1:C-1,2)));
 356
 357
 358 if length (col) == 1            # Convert graylevel to RGB
 359   col = [1 1 1]*col;
 360
 361 elseif any (prod (size (col)) == [R*C,(R-1)*(C-1)])
 362   col = [1;1;1]*col(:)';
 363 end
 364
 365 if zgray || zrb || any (zcol(:)) # Treat zgray zrb and zcol options
 366   zx = max (z(keepip));
 367   zn = min (z(keepip));
 368   if     zgray, zcol = [0 0 0; 1 1 1]';
 369   elseif zrb  , zcol = [0 0 0.7; 0.5 0 0.8; 1 0 0]';
 370   end
 371
 372   ci = 1 + floor (cw = (columns (zcol)-1) * (z(keepip) - zn)/(zx - zn));
 373   cw =  cw - ci + 1;
 374   ii = find (ci >= columns (zcol));
 375   if ! isempty (ii), ci(ii) = columns (zcol) - 1; cw(ii) = 1; end
 376   col = zeros (3,R*C);
 377   col(:,keepip) = \
 378       zcol(:,ci) .* ([1;1;1]*(1-cw)) + zcol(:,ci+1) .* ([1;1;1]*cw);
 379
 380 end                             # EOF zgray zrb and zcol options
 381
 382
 383 if checker && numel (col) <= 6
 384   if isnan (checker), checker = 10; end
 385   if length (checker) == 1, checker = [checker, checker]; end
 386
 387   col = checker_color (checker, col, R,C);
 388
 389   if 0
 390     if checker(1) > 0, checker(1) = - (C-1)/checker(1); end
 391     if checker(2) > 0, checker(2) = - (R-1)/checker(2); end
 392
 393     checker *= -1;
 394     colx = 2 * (rem (0:C-2,2*checker(1)) < checker(1)) - 1;
 395     coly = 2 * (rem (0:R-2,2*checker(2)) < checker(2)) - 1;
 396     icol = 1 + ((coly'*colx) > 0);
 397                                 # Keep at most 1st 2 colors of col for the
 398                                 # checker
 399     if prod (size (col)) == 2,
 400       col = [1;1;1]*col;
 401     elseif  prod (size (col)) < 6, # Can't be < 3 because of previous code
 402       col = col(1:3)(:);
 403       if all (col >= 1-eps), col = [col [0;0;0]];       # Black and White
 404       else                   col = [col [1;1;1]];       # X and White
 405       end
 406     end
 407     col = reshape (col(:),3,2);
 408     col = col(:,icol);
 409   end
 410 end                             # EOF if checker
 411
 412
 413 if prod (size (col)) == 3*(R-1)*(C-1),
 414   colorPerVertex = 0;
 415 end
 416
 417 if ! colorPerVertex
 418   if prod (size (col)) == 3*(R-1)*(C-1)
 419     col = reshape (col,3, (R-1)*(C-1));
 420     col = [col, col];
 421   else
 422     printf(["vrml_surf : ",\
 423             " colorPerVertex==0, (R-1)*(C-1)==%i, but col has size [%i,%i]\n"],\
 424      R*C,size (col));
 425   end
 426
 427 end
 428
 429 if ! all(keepp),
 430
 431                                 # Try to toggle some triangles to fill in
 432                                 # holes
 433   if ! upper
 434     nt = (R-1)*(C-1) ;
 435     tmp = ! reshape(keepp(trgs),3,2*nt);
 436     tmp = all( tmp == kron([0,0;0,1;1,0],ones(1,nt)) );
 437     trgs(3,     find (    tmp(1:nt)   & rem (trgs(3,1:nt),R)) )++ ;
 438     trgs(2, nt+ find ( tmp(nt+1:2*nt) & (rem (trgs(3,nt+1:2*nt),R)!=1)) )-- ;
 439
 440                                 # Remove whatever can't be kept
 441     keept = all (reshape(keepp(trgs),3,2*(R-1)*(C-1)));
 442   else
 443     tmp = reshape (keepp,R,C);
 444     keept = \
 445         all (reshape (tmp(trgs(1:2,:)),2,2*(R-1)*(C-1))) & \
 446         [(tmp(1:R-1,2:C) | tmp(2:R,2:C))(:)', \
 447          (tmp(1:R-1,1:C-1) | tmp(2:R,1:C-1))(:)'] ;
 448   end
 449
 450   keepit = find (keept);
 451
 452   renum = cumsum (keepp);
 453
 454   pts = pts (:,keepip);
 455   trgs = reshape(renum (trgs (:,keepit)), 3, columns(keepit));
 456
 457   if prod (size (col)) == 6*(R-1)*(C-1)
 458     col = col(:,keepit);
 459                                 # Coherence check : colorPerVertex == 0
 460     if colorPerVertex
 461       error ("Col has size 3*(R-1)*(C-1), but colorPerVertex == 1");
 462     end
 463
 464   elseif prod (size (col)) == 3*R*C
 465     col = col(:,keepip);
 466
 467                                 # Coherence check : colorPerVertex == 1
 468     if ! colorPerVertex
 469       error ("Col has size 3*R*C, but colorPerVertex == 0");
 470     end
 471   end
 472
 473 end
 474 ## printf ("Calling vrml_faces\n");
 475 if !tex
 476   s = vrml_faces (pts, trgs, "col", col,\
 477                   "colorPerVertex",colorPerVertex,\
 478                   "creaseAngle", creaseAngle,\
 479                   "tran", tran, "emit", emit,\
 480                   "DEFcoord",DEFcoord,"DEFcol",DEFcol);
 481 else
 482    texXY = texXY(:,keepip);
 483 #   texXY(:,[1:5,232:236])
 484 #   pts(:,[1:5,232:236])
 485 #   trgs(:,1:20)
 486 #  [texXY;  pts]
 487 #  trgs
 488 #  texXY(:,trgs(:))
 489 #   R, C
 490 #  keyboard
 491   s = vrml_faces (pts, trgs,\
 492                   "tran", tran, "tex", tex, "tcoord", texXY,\
 493                   "DEFcoord",DEFcoord,"DEFcol",DEFcol);
 494 end
 495 ## printf ("Done\n");
 496 ## R=5; C=11;
 497 ## x = ones(R,1)*[1:C]; y = [1:R]'*ones(1,C);
 498 ## zz = z = sin(x)+(2*y/R-1).^2 ;
 499 ## ## Punch some holes
 500 ## holes = ind2sub ([R,C],[3,3,3,1,2,3,2,3;1,2,3,5,7,7,9,10]')
 501 ## z(holes) = nan
 502 ## save_vrml("tmp.wrl",vrml_surf(x,y,z+1))
 503 ## save_vrml("tmp.wrl",vrml_surf(z,"col",[0.5,0,0],"tran",0.5),vrml_surf(z+1))
 504
 505 endfunction
 506