POV-Ray : Newsgroups : povray.text.scene-files : Code for NURBS surfaces as mesh2{} objects : Code for NURBS surfaces as mesh2{} objects Server Time
25 Apr 2024 11:21:50 EDT (-0400)
  Code for NURBS surfaces as mesh2{} objects  
From: Tor Olav Kristensen
Date: 25 Mar 2002 15:28:38
Message: <3C9F86D5.778C126F@hotmail.com>
A little bit later than promised, but here it is.

I recommend strongly to read a little about NURBS
and POV functions before using the code below.

Note that the code below will not work with POV-
Ray v3.5 Beta 13.

This is work in progress, so there might be many
errors and that anything may change.

I'll explain more about how to use the macros later
(if I can find time).

Now have fun.


Tor Olav


// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Copyright 2002 by Tor Olav Kristensen
// Email: tor### [at] hotmailcom
// http://www.crosswinds.net/~tok
// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

// A simple NURBS surface as a mesh2{} object

#version 3.5; // Beta 14 or later (?)

#include "colors.inc"

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Macros needed because of bug in POV-Ray v3.1 and v3.5 Beta 1 to 13

#macro writeM(AA, FileName)

  #local Dim = dimension_size(AA, 1);
  #fopen ArrayFile FileName write
  #write (ArrayFile, "array[", Dim, "] { ")
  #local Cnt = 0;
  #while (Cnt < Dim)
    #write (ArrayFile, AA[Cnt])
    #if (Cnt < Dim - 1)
      #write (ArrayFile, ", ")
    #end // if
    #local Cnt = Cnt + 1;
  #end // while
  #write (ArrayFile, " }\n")
  #fclose ArrayFile

#end // macro writeM


#macro Mread(FileName)

  #include FileName

#end // macro Mread


#macro Mret(AA)

  writeM(AA, "Temp.txt")
  Mread("Temp.txt")

#end // macro Mret


#macro writeMM(AA, FileName)

  #local DimR = dimension_size(AA, 1);
  #local DimC = dimension_size(AA, 2);
  #fopen ArrayFile FileName write
  #write (ArrayFile, "array[", DimR, "][", DimC,"] {\n")
  #local CntR = 0;
  #while (CntR < DimR)
    #write (ArrayFile, " { ")
    #local CntC = 0;
    #while (CntC < DimC)
      #write (ArrayFile, AA[CntR][CntC])
      #if (CntC < DimC - 1)
        #write (ArrayFile, ", ")
      #end // if
      #local CntC = CntC + 1;
    #end // while
    #write (ArrayFile, " }")
    #if (CntR < DimR - 1)
      #write (ArrayFile, ",")
    #end // if
    #write (ArrayFile, "\n")
    #local CntR = CntR + 1;
  #end // while
  #write (ArrayFile, "}\n")
  #fclose ArrayFile

#end // macro writeMM


#macro MMread(FileName)

  #include FileName

#end // macro MMread


#macro MMret(AA)

  writeMM(AA, "Temp.txt")
  MMread("Temp.txt")

#end // macro MMret

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

#macro printM(AA)

  #local Dim = dimension_size(AA, 1);
  #debug concat("array[", str(Dim, 0, 0), "] { ")
  #local Cnt = 0;
  #while (Cnt < Dim)
    #debug str(AA[Cnt], 0, -1)
    #if (Cnt < Dim - 1)
      #debug ", "
    #end // if
    #local Cnt = Cnt + 1;
  #end // while
  #debug " }\n"

#end // macro printMM


#macro printMM(AA)

  #local DimR = dimension_size(AA, 1);
  #local DimC = dimension_size(AA, 2);
  #debug concat("array[", str(DimR, 0, 0), "]")
  #debug concat("[", str(DimC,0, 0), "] {\n")
  #local CntR = 0;
  #while (CntR < DimR)
    #debug " { "
    #local CntC = 0;
    #while (CntC < DimC)
      #debug str(AA[CntR][CntC], 0, -1)
      #if (CntC < DimC - 1)
        #debug ", "
      #end // if
      #local CntC = CntC + 1;
    #end // while
    #debug " }"
    #if (CntR < DimR - 1)
      #debug ","
    #end // if
    #debug "\n"
    #local CntR = CntR + 1;
  #end // while
  #debug "}\n"

#end // macro printMM

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Macros for making a mesh from 3 uv-functions

#macro ParametricMesh(xxFn, yyFn, zzFn, UVmin, UVmax, NrUV)

  #local NrU = NrUV.u;
  #local NrV = NrUV.v;
  #local NrUP = NrU + 1;
  #local NrVP = NrV + 1;
  #local NrOfVertices = NrUP*NrVP;
  #local Nothing = 1E-9;
  #local dUV = (UVmax - UVmin - Nothing*<1, 1>)/<NrU, NrV>; 

  #local pVertices = array[NrUP][NrVP]

  #local UU = 0;
  #while (UU < NrUP)
    #local VV = 0;
    #while (VV < NrVP)
      #local UV = UVmin + <UU, VV>*dUV;
      #local U = UV.u;
      #local V = UV.v;
      #local pVertices[UU][VV] = <xxFn(U, V), yyFn(U, V), zzFn(U, V)>;
      #local VV = VV + 1;
    #end // while
    #local UU = UU + 1;
  #end // while

  #local NrOfFaces =2*NrU*NrV;
  #local vFaces1 = array[NrU][NrV]
  #local vFaces2 = array[NrU][NrV]
  #local vFace1Normals = array[NrU][NrV]
  #local vFace2Normals = array[NrU][NrV]
  #local FST = array[NrU][NrV] // Face Split Type

  #local UU = 0;
  #while (UU < NrU)
    #local UP = UU + 1;
    #local VV = 0;
    #while (VV < NrV)
      #local VP = VV + 1;
      #local UVmid = UVmin + <UU + 0.5, VV + 0.5>*dUV;
      #local U = UVmid.u;
      #local V = UVmid.v;
      #local pMid = <xxFn(U, V), yyFn(U, V), zzFn(U, V)>;
      #local pUUVV = pVertices[UU][VV];
      #local pUPVP = pVertices[UP][VP];
      #local pUPVV = pVertices[UP][VV];
      #local pUUVP = pVertices[UU][VP];
      #local pMeanF = (pUUVV + pUPVP)/2;
      #local pMeanB = (pUPVV + pUUVP)/2;
      #if (vlength(pMeanF - pMid) < vlength(pMeanB - pMid))
        #local FST[UU][VV] = 0;
        #local vFace1Normals[UU][VV] =
          vcross(pUPVP - pUUVV, pUUVP - pUUVV);
        #local vFace2Normals[UU][VV] =
          vcross(pUPVV - pUUVV, pUPVP - pUUVV);
        #local vFaces1[UU][VV] = <UU, UP, UU>*NrVP + <VV, VP, VP>;
        #local vFaces2[UU][VV] = <UU, UP, UP>*NrVP + <VV, VV, VP>;
      #else
        #local FST[UU][VV] = 1;
        #local vFace1Normals[UU][VV] =
          vcross(pUPVV - pUUVV, pUUVP - pUUVV);
        #local vFace2Normals[UU][VV] =
          vcross(pUUVP - pUPVP, pUPVV - pUPVP);
        #local vFaces1[UU][VV] = <UU, UP, UU>*NrVP + <VV, VV, VP>;
        #local vFaces2[UU][VV] = <UP, UU, UP>*NrVP + <VP, VP, VV>;
      #end // if
      #local VV = VP;
    #end // while
    #local UU = UP;
  #end // while

  #local vNormals = array[NrUP][NrVP]

  #local UU = 0;
  #while (UU < NrUP)
    #local NFU = (UU > 0); // Not First U
    #local NLU = (UU < NrU); // Not Last U
    #local UM = UU - 1; // U Minus 1
    #local VV = 0;
    #while (VV < NrVP)
      #local NFV = (VV > 0); // Not First V
      #local NLV = (VV < NrV); // Not Last V
      #local VM = VV - 1; // V Minus 1
      #local vNormals[UU][VV] =
        vnormalize(
          <0, 0, 0>
          #if (NFU)
            #if (NFV)
              #if (FST[UM][VM] = 0)
               +vFace1Normals[UM][VM]
              #end // if
             +vFace2Normals[UM][VM]
            #end // if
            #if (NLV)
              #if (FST[UM][VV] = 1)
               +vFace1Normals[UM][VV]
              #end // if
             +vFace2Normals[UM][VV]
            #end // if
          #end // if
          #if (NLU)
            #if (NFV)
             +vFace1Normals[UU][VM]
              #if (FST[UU][VM] = 1)
               +vFace2Normals[UU][VM]
              #end // if
            #end // if
            #if (NLV)
             +vFace1Normals[UU][VV]
              #if (FST[UU][VV] = 0)
               +vFace2Normals[UU][VV]
              #end // if
            #end // if
          #end // if
        );
      #local VV = VV + 1;
    #end // while
    #local UU = UU + 1;
  #end // while

  mesh2 {
    vertex_vectors {
      NrOfVertices,
      #local UU = 0;
      #while (UU < NrUP)
        #local VV = 0;
        #while (VV < NrVP)
          pVertices[UU][VV],
          #local VV = VV + 1;
        #end // while
        #local UU = UU + 1;
      #end // while
    }
    uv_vectors {
      NrOfVertices,
      #local UU = 0;
      #while (UU < NrUP)
        #local VV = 0;
        #while (VV < NrVP)
          (UVmin + <UU, VV>*dUV),
          #local VV = VV + 1;
        #end // while
        #local UU = UU + 1;
      #end // while
    }
    normal_vectors {
      NrOfVertices,
      #local UU = 0;
      #while (UU < NrUP)
        #local VV = 0;
        #while (VV < NrVP)
          vNormals[UU][VV],
          #local VV = VV + 1;
        #end // while
        #local UU = UU + 1;
      #end // while
    }
    face_indices {
      NrOfFaces,
      #local UU = 0;
      #while (UU < NrU)
        #local VV = 0;
        #while (VV < NrV)
          vFaces1[UU][VV],
          vFaces2[UU][VV],
          #local VV = VV + 1;
        #end // while
        #local UU = UU + 1;
      #end // while
    }
  }

#end // macro ParametricMesh


#macro WriteParamMesh(xxFn, yyFn, zzFn, UVmin, UVmax, NrUV, FileName)

  #local NrU = NrUV.u;
  #local NrV = NrUV.v;
  #local NrUP = NrU + 1;
  #local NrVP = NrV + 1;
  #local NrOfVertices = NrUP*NrVP;
  #local Nothing = 1E-9;
  #local dUV = (UVmax - UVmin - Nothing*<1, 1>)/<NrU, NrV>; 

  #local pVertices = array[NrUP][NrVP]

  #local UU = 0;
  #while (UU < NrUP)
    #local VV = 0;
    #while (VV < NrVP)
      #local UV = UVmin + <UU, VV>*dUV;
      #local U = UV.u;
      #local V = UV.v;
      #local pVertices[UU][VV] = <xxFn(U, V), yyFn(U, V), zzFn(U, V)>;
      #local VV = VV + 1;
    #end // while
    #local UU = UU + 1;
  #end // while

  #local NrOfFaces =2*NrU*NrV;
  #local vFaces1 = array[NrU][NrV]
  #local vFaces2 = array[NrU][NrV]
  #local vFace1Normals = array[NrU][NrV]
  #local vFace2Normals = array[NrU][NrV]
  #local FST = array[NrU][NrV] // Face Split Type

  #local UU = 0;
  #while (UU < NrU)
    #local UP = UU + 1;
    #local VV = 0;
    #while (VV < NrV)
      #local VP = VV + 1;
      #local UVmid = UVmin + <UU + 0.5, VV + 0.5>*dUV;
      #local U = UVmid.u;
      #local V = UVmid.v;
      #local pMid = <xxFn(U, V), yyFn(U, V), zzFn(U, V)>;
      #local pUUVV = pVertices[UU][VV];
      #local pUPVP = pVertices[UP][VP];
      #local pUPVV = pVertices[UP][VV];
      #local pUUVP = pVertices[UU][VP];
      #local pMeanF = (pUUVV + pUPVP)/2;
      #local pMeanB = (pUPVV + pUUVP)/2;
      #if (vlength(pMeanF - pMid) < vlength(pMeanB - pMid))
        #local FST[UU][VV] = 0;
        #local vFace1Normals[UU][VV] =
          vcross(pUPVP - pUUVV, pUUVP - pUUVV);
        #local vFace2Normals[UU][VV] =
          vcross(pUPVV - pUUVV, pUPVP - pUUVV);
        #local vFaces1[UU][VV] = <UU, UP, UU>*NrVP + <VV, VP, VP>;
        #local vFaces2[UU][VV] = <UU, UP, UP>*NrVP + <VV, VV, VP>;
      #else
        #local FST[UU][VV] = 1;
        #local vFace1Normals[UU][VV] =
          vcross(pUPVV - pUUVV, pUUVP - pUUVV);
        #local vFace2Normals[UU][VV] =
          vcross(pUUVP - pUPVP, pUPVV - pUPVP);
        #local vFaces1[UU][VV] = <UU, UP, UU>*NrVP + <VV, VV, VP>;
        #local vFaces2[UU][VV] = <UP, UU, UP>*NrVP + <VP, VP, VV>;
      #end // if
      #local VV = VP;
    #end // while
    #local UU = UP;
  #end // while

  #local vNormals = array[NrUP][NrVP]

  #local UU = 0;
  #while (UU < NrUP)
    #local NFU = (UU > 0); // Not First U
    #local NLU = (UU < NrU); // Not Last U
    #local UM = UU - 1; // U Minus 1
    #local VV = 0;
    #while (VV < NrVP)
      #local NFV = (VV > 0); // Not First V
      #local NLV = (VV < NrV); // Not Last V
      #local VM = VV - 1; // V Minus 1
      #local vNormals[UU][VV] =
        vnormalize(
          <0, 0, 0>
          #if (NFU)
            #if (NFV)
              #if (FST[UM][VM] = 0)
               +vFace1Normals[UM][VM]
              #end // if
             +vFace2Normals[UM][VM]
            #end // if
            #if (NLV)
              #if (FST[UM][VV] = 1)
               +vFace1Normals[UM][VV]
              #end // if
             +vFace2Normals[UM][VV]
            #end // if
          #end // if
          #if (NLU)
            #if (NFV)
             +vFace1Normals[UU][VM]
              #if (FST[UU][VM] = 1)
               +vFace2Normals[UU][VM]
              #end // if
            #end // if
            #if (NLV)
             +vFace1Normals[UU][VV]
              #if (FST[UU][VV] = 0)
               +vFace2Normals[UU][VV]
              #end // if
            #end // if
          #end // if
        );
      #local VV = VV + 1;
    #end // while
    #local UU = UU + 1;
  #end // while

  #fopen Meshfile FileName write
  #write (Meshfile, "mesh2 {\n")
  #write (Meshfile, "  vertex_vectors {\n")
  #write (Meshfile, "    ", NrOfVertices, ",\n")
  #local UU = 0;
  #while (UU < NrUP)
    #local VV = 0;
    #while (VV < NrVP)
      #write (Meshfile, "    <")
      #write (Meshfile, vstr(3, pVertices[UU][VV], ", ", 0, -1))
      #write (Meshfile, ">,\n")
      #local VV = VV + 1;
    #end // while
    #local UU = UU + 1;
  #end // while
  #write (Meshfile, "  }\n")
  #write (Meshfile, "  uv_vectors {\n")
  #write (Meshfile, "    ", NrOfVertices, ",\n")
  #local UU = 0;
  #while (UU < NrUP)
    #local VV = 0;
    #while (VV < NrVP)
      #write (Meshfile, "    <")
      #write (Meshfile, vstr(2, UVmin + <UU, VV>*dUV, ", ", 0, -1))
      #write (Meshfile, ">,\n")
      #local VV = VV + 1;
    #end // while
    #local UU = UU + 1;
  #end // while
  #write (Meshfile, "  }\n")
  #write (Meshfile, "  normal_vectors {\n")
  #write (Meshfile, "    ", NrOfVertices, ",\n")
  #local UU = 0;
  #while (UU < NrUP)
    #local VV = 0;
    #while (VV < NrVP)
      #write (Meshfile, "    <")
      #write (Meshfile, vstr(3, vNormals[UU][VV], ", ", 0, -1))
      #write (Meshfile, ">,\n")
      #local VV = VV + 1;
    #end // while
    #local UU = UU + 1;
  #end // while
  #write (Meshfile, "  }\n")
  #write (Meshfile, "  face_indices {\n")
  #write (Meshfile, "    ", NrOfFaces, ",\n")
  #local UU = 0;
  #while (UU < NrU)
    #local VV = 0;
    #while (VV < NrV)
      #write (Meshfile, "    <")
      #write (Meshfile, vstr(3, vFaces1[UU][VV], ", ", 0, -1))
      #write (Meshfile, ">, <")
      #write (Meshfile, vstr(3, vFaces2[UU][VV], ", ", 0, -1))
      #write (Meshfile, ">,\n")
      #local VV = VV + 1;
    #end // while
    #local UU = UU + 1;
  #end // while
  #write (Meshfile, "  }\n")
  #write (Meshfile, "}\n")
  #fclose Meshfile

#end // macro WriteParamMesh


#macro CylinderMesh(xxFn, yyFn, zzFn, UVmin, UVmax, NrUV, CylRad)

  #local NrU = NrUV.u;
  #local NrV = NrUV.v;
  #local NrUP = NrU + 1;
  #local NrVP = NrV + 1;
  #local NrOfVertices = NrUP*NrVP;
  #local Nothing = 1E-9;
  #local dUV = (UVmax - UVmin - Nothing*<1, 1>)/<NrU, NrV>; 

  #local pVertices = array[NrUP][NrVP]
  #local UU = 0;
  #while (UU < NrUP)
    #local VV = 0;
    #while (VV < NrVP)
      #local UV = UVmin + <UU, VV>*dUV;
      #local U = UV.u;
      #local V = UV.v;
      #local pVertices[UU][VV] = <xxFn(U, V), yyFn(U, V), zzFn(U, V)>;
      #local VV = VV + 1;
    #end // while
    #local UU = UU + 1;
  #end // while

  union {
    #local UU = 0;
    #while (UU < NrU)
      #local UP = UU + 1;
      #local VV = 0;
      #while (VV < NrV)
        #local VP = VV + 1;
        #local UVmid = UVmin + <UU + 0.5, VV + 0.5>*dUV;
        #local U = UVmid.u;
        #local V = UVmid.v;
        #local pMid = <xxFn(U, V), yyFn(U, V), zzFn(U, V)>;
        #local pUUVV = pVertices[UU][VV];
        #local pUPVP = pVertices[UP][VP];
        #local pUPVV = pVertices[UP][VV];
        #local pUUVP = pVertices[UU][VP];
        #local pMeanF = (pUUVV + pUPVP)/2;
        #local pMeanB = (pUPVV + pUUVP)/2;
        #local LF = vlength(pMeanF - pMid);
        #local LB = vlength(pMeanB - pMid);
        #if (LF < LB)
          cylinder {
            pUUVV, pUPVP, CylRad
//            pigment { color Red }
          }
        #else
          cylinder {
            pUPVV, pUUVP, CylRad
//            pigment { color Red }
          }
        #end // if
        cylinder {
          pUUVV, pUPVV, CylRad
//          pigment { color White }
        }
        cylinder {
          pUUVV, pUUVP, CylRad
//          pigment { color White }
        }
        #if (VV = NrV - 1)
          cylinder {
            pUUVP, pUPVP, CylRad
//            pigment { color White }
          }
        #end // if
        #if (UU = NrU - 1)
          cylinder {
            pUPVV, pUPVP, CylRad
//            pigment { color White }
          }
        #end // if
        #local VV = VP;
      #end // while
      #local UU = UP;
    #end // while
  }

#end // macro CylinderMesh

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Macros for building blending functions
// v3.5 Beta 14 (or later ?) is needed for these macros

#macro DeCasteljau_U(Order, Knots, K)

  #local tA = Knots[K];
  #local tB = Knots[K + 1];
  #local tC = Knots[K + Order];
  #local tD = Knots[K + Order - 1];
  #local d1 = tD - tA;
  #local d2 = tC - tB;
  
  #if (Order > 1)
    #if (d1 = 0)
      #if (d2 = 0)
        0
      #else
        #if (tC = 0)
          -u
        #else
          (tC - u)
        #end // if
        #if (d2 != 1)
          /d2
        #end // if
        *DeCasteljau_U(Order - 1, Knots, K + 1)
      #end // if
    #else
      #if (d2 = 0)
        #if (tA = 0)
          u
        #else
          (u - tA)
        #end // if
        #if (d1 != 1)
          /d1
        #end // if
        *DeCasteljau_U(Order - 1, Knots, K)
      #else
        (
          #if (tA = 0)
            u
          #else
            (u - tA)
          #end // if
          #if (d1 != 1)
            /d1
          #end // if
          *DeCasteljau_U(Order - 1, Knots, K)
          #if (tC = 0)
            -u
          #else
            +(tC - u)
          #end // if
          #if (d2 != 1)
            /d2
          #end // if
          *DeCasteljau_U(Order - 1, Knots, K + 1)
        )
      #end // if
    #end // if
  #else
    #if (tA = 0)
      #if (tB = 0)
        select(u, 0, 1, 0) // or 0 ?
      #else
        select(u, 0, select(u - tB, 1, 0))
      #end // if
    #else
      #if (tB = 0)
        select(u - tA, 0, select(u, 1, 0))
      #else
        #if (tA = tB)
          select(u - tA, 0, 1, 0) // or 0 ?
        #else
          select(u - tA, 0, select(u - tB, 1, 0))
        #end // if
      #end // if
    #end // if
  #end // if

#end // macro DeCasteljau_U


#macro DeCasteljau_V(Order, Knots, K)

  #local tA = Knots[K];
  #local tB = Knots[K + 1];
  #local tC = Knots[K + Order];
  #local tD = Knots[K + Order - 1];
  #local d1 = tD - tA;
  #local d2 = tC - tB;
  
  #if (Order > 1)
    #if (d1 = 0)
      #if (d2 = 0)
        0
      #else
        #if (tC = 0)
          -v
        #else
          (tC - v)
        #end // if
        #if (d2 != 1)
          /d2
        #end // if
        *DeCasteljau_V(Order - 1, Knots, K + 1)
      #end // if
    #else
      #if (d2 = 0)
        #if (tA = 0)
          v
        #else
          (v - tA)
        #end // if
        #if (d1 != 1)
          /d1
        #end // if
        *DeCasteljau_V(Order - 1, Knots, K)
      #else
        (
          #if (tA = 0)
            v
          #else
            (v - tA)
          #end // if
          #if (d1 != 1)
            /d1
          #end // if
          *DeCasteljau_V(Order - 1, Knots, K)
          #if (tC = 0)
            -v
          #else
            +(tC - v)
          #end // if
          #if (d2 != 1)
            /d2
          #end // if
          *DeCasteljau_V(Order - 1, Knots, K + 1)
        )
      #end // if
    #end // if
  #else
    #if (tA = 0)
      #if (tB = 0)
        select(v, 0, 1, 0) // or 0 ?
      #else
        select(v, 0, select(v - tB, 1, 0))
      #end // if
    #else
      #if (tB = 0)
        select(v - tA, 0, select(v, 1, 0))
      #else
        #if (tA = tB)
          select(v - tA, 0, 1, 0) // or 0 ?
        #else
          select(v - tA, 0, select(v - tB, 1, 0))
        #end // if
      #end // if
    #end // if
  #end // if

#end // macro DeCasteljau_V


#macro BlendUFunction(Order, Knots, K)

  function(u) { DeCasteljau_U(Order, Knots, K) }

#end // macro BlendUFunction


#macro BlendVFunction(Order, Knots, K)

  function(v) { DeCasteljau_V(Order, Knots, K) }

#end // macro BlendVFunction

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Macros below here are useable for "direct" user calls

#macro NURBS_SurfFunction(Order, KnotsU, KnotsV, Points, Weights, v0)

  #local Usize = dimension_size(Points, 1);
  #local Vsize = dimension_size(Points, 2);
  #local BlendUFn = array[Usize]
  #local BlendVFn = array[Vsize]
  #local I = 0;
  #while (I < Usize)
    #local BlendUFn[I] = BlendUFunction(Order, KnotsU, I)
    #local I = I + 1;
  #end // while
  #local J = 0;
  #while (J < Vsize)
    #local BlendVFn[J] = BlendVFunction(Order, KnotsV, J)
    #local J = J + 1;
  #end // while

  #local TFn =
    function(u, v) {
      #local I = 0;
      #while (I < Usize)
        +BlendUFn[I](u)*
        (
          0
          #local J = 0;
          #while (J < Vsize)
            #local WC = Weights[I][J]*vdot(Points[I][J], v0);
            #if (WC != 0)
              +BlendVFn[J](v)
              #if (WC != 1)
                *WC
              #end // if
            #end // if
            #local J = J + 1;
          #end // while
        )
        #local I = I + 1;
      #end // while
    }

  #local NFn =
    function(u, v) {
      #local I = 0;
      #while (I < Usize)
        +BlendUFn[I](u)*
        (
          0
          #local J = 0;
          #while (J < Vsize)
            #local W = Weights[I][J];
            #if (W != 0)
              +BlendVFn[J](v)
              #if (W != 1)
                *W
              #end // if
            #end // if
            #local J = J + 1;
          #end // while
        )
        #local I = I + 1;
      #end // while
    }

  function(u, v) { TFn(u, v)/NFn(u, v) }

#end // macro NURBS_SurfFunction


#macro OpenUniformKnotVector(Order, Nr, Tmin, Tmax)

  #local Size = Order + Nr;
  #local KK = array[Size]
  #local dT = (Tmax - Tmin)/(Size - 2*Order + 1);
  #local Cnt = 0;
  #while (Cnt < Size)
    #local T = (Cnt - Order + 1)*dT;
    #local KK[Cnt] = min(max(Tmin, T), Tmax);
    #local Cnt = Cnt + 1;
  #end // while

  Mret(KK)

#end // macro OpenUniformKnotVector


#macro UnitWeights(SizeU, SizeV)

  #local Weights = array[SizeU][SizeV]
  #local CntU = 0;
  #while (CntU < SizeU)
    #local CntV = 0;
    #while (CntV < SizeV)
      #local Weights[CntU][CntV] = 1;
      #local CntV = CntV + 1;
    #end // while
    #local CntU = CntU + 1;
  #end // while

  MMret(Weights)

#end // macro UnitWeights


#macro NURBS_SurfaceMesh(Points, MeshSize)

  #local Order = 4;
  #local SizeU = dimension_size(Points, 1);
  #local SizeV = dimension_size(Points, 2);
  #local MinUV = <0, 0>;
  #local MaxUV = <1, 1>;
  #local KnotsU = OpenUniformKnotVector(Order, SizeU, MinUV.u, MaxUV.u)
  #local KnotsV = OpenUniformKnotVector(Order, SizeV, MinUV.v, MaxUV.v)
  #local Weights = UnitWeights(SizeU, SizeV)

  ParametricMesh(
    NURBS_SurfFunction(Order, KnotsU, KnotsV, Points, Weights, x),
    NURBS_SurfFunction(Order, KnotsU, KnotsV, Points, Weights, y),
    NURBS_SurfFunction(Order, KnotsU, KnotsV, Points, Weights, z),
    MinUV,
    MaxUV,
    MeshSize
  )

#end // NURBS_SurfaceMesh


#macro NURBS_CylinderMesh(Points, MeshSize, CylinderRadius)

  #local Order = 4;
  #local SizeU = dimension_size(Points, 1);
  #local SizeV = dimension_size(Points, 2);
  #local MinUV = <0, 0>;
  #local MaxUV = <1, 1>;
  #local KnotsU = OpenUniformKnotVector(Order, SizeU, MinUV.u, MaxUV.u)
  #local KnotsV = OpenUniformKnotVector(Order, SizeV, MinUV.v, MaxUV.v)
  #local Weights = UnitWeights(SizeU, SizeV)

  CylinderMesh(
    NURBS_SurfFunction(Order, KnotsU, KnotsV, Points, Weights, x),
    NURBS_SurfFunction(Order, KnotsU, KnotsV, Points, Weights, y),
    NURBS_SurfFunction(Order, KnotsU, KnotsV, Points, Weights, z),
    MinUV,
    MaxUV,
    MeshSize,
    CylinderRadius
  )

#end // NURBS_CylinderMesh


#macro ShowControlPoints(Points, Radius)

  union {
    #local SizeU = dimension_size(Points, 1);
    #local SizeV = dimension_size(Points, 2);
    #local CntU = 0;
    #while (CntU < SizeU)
      #local CntV = 0;
      #while (CntV < SizeV)
        sphere { Points[CntU][CntV], Radius }
        #local CntV = CntV + 1;
      #end // while
      #local CntU = CntU + 1;
    #end // while
  }

#end // macro ShowControlPoints


#macro ShowControlGrid(Points, Radius)

  #local SizeU = dimension_size(Points, 1);
  #local SizeV = dimension_size(Points, 2);

  union {
    #local CntU = 0;
    #while (CntU < SizeU - 1)
      #local CntV = 0;
      #while (CntV < SizeV - 1)
        cylinder {
          Points[CntU][CntV],
          Points[CntU + 1][CntV],
          Radius
        }
        cylinder {
          Points[CntU][CntV],
          Points[CntU][CntV + 1],
          Radius
        }
        #local CntV = CntV + 1;
      #end // while
      #local CntU = CntU + 1;
    #end // while
  
    #local CntU = 0;
    #while (CntU < SizeU - 1)
      cylinder {
        Points[CntU][SizeV - 1],
        Points[CntU + 1][SizeV - 1],
        Radius
      }
      #local CntU = CntU + 1;
    #end // while
  
    #local CntV = 0;
    #while (CntV < SizeV - 1)
      cylinder {
        Points[SizeU - 1][CntV],
        Points[SizeU - 1][CntV + 1],
        Radius
      }
      #local CntV = CntV + 1;
    #end // while
  }

#end // macro ShowControlGrid

// ===== 1 ================= 3 ================= 5 ================= 7
// =============== 2 ================= 4 ================= 6 =========
// Now lets make a NURBS surface . . .

// Define some points to control the NURBS surface

#declare NrU = 9;
#declare NrV = 7;

#declare CtrlPoints =
  array[NrU][NrV] {
    { < 6,-9,-9>, < 6,-9,-8>, < 6,-9,-7>, < 6,-9,-6>, 
                              < 7,-9,-6>, < 8,-9,-6>, < 9,-9,-6> },
    { < 3,-9,-9>, < 3,-9,-7>, < 3,-9,-5>, < 3,-9,-3>,
                              < 5,-9,-3>, < 7,-9,-3>, < 9,-9,-3> },
    { < 0,-9,-9>, < 0,-9,-6>, < 0,-9,-3>, < 0,-9, 0>,
                              < 3,-9, 0>, < 6,-9, 0>, < 9,-9, 0> },
    { < 0,-6,-9>, < 0,-6,-6>, < 0,-6,-3>, < 0,-6, 0>,
                              < 3,-6, 0>, < 6,-6, 0>, < 9,-6, 0> },
    { < 0,-3,-9>, < 0,-3,-6>, < 0,-3,-3>, < 0,-3, 0>,
                              < 3,-3, 0>, < 6,-3, 0>, < 9,-3, 0> },
    { < 0, 0,-9>, < 0, 0,-6>, < 0, 0,-3>, < 0, 0, 0>,
                              < 3, 0, 0>, < 6, 0, 0>, < 9, 0, 0> },
    { <-3, 0,-9>, <-3, 0,-5>, <-3, 0,-1>, <-3, 0, 3>,
                              < 1, 0, 3>, < 5, 0, 3>, < 9, 0, 3> },
    { <-6, 0,-9>, <-6, 0,-4>, <-6, 0, 1>, <-6, 0, 6>,
                              <-1, 0, 6>, < 4, 0, 6>, < 9, 0, 6> },
    { <-9, 0,-9>, <-9, 0,-3>, <-9, 0, 3>, <-9, 0, 9>,
                              <-3, 0, 9>, < 3, 0, 9>, < 9, 0, 9> }
  }

/*
object {
  ShowControlPoints(CtrlPoints, 0.16)
  pigment { color Yellow }
}

object {
  ShowControlGrid(CtrlPoints, 0.08)
  pigment { color Cyan }
}
*/

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// And then call two macros

object {
  NURBS_SurfaceMesh(CtrlPoints, <50, 50>)
  pigment { color Blue + Cyan/2 }
}

object {
  NURBS_CylinderMesh(CtrlPoints, <50, 50>, 0.02)
  pigment { color Magenta/2 + White/4 + Red/2 }
}

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// This does the same, but should parse in a little less time

/*
#declare Order = 4;

#declare Weights = UnitWeights(NrU, NrV)

// Uncomment the 2 lines below to see the surface deformed locally
//#declare CtrlPoints[4][2] = CtrlPoints[4][2] + <12, 12, -6>;
//#declare Weights[4][2] = 0.5;

#declare MinUV = <0, 0>;
#declare MaxUV = <1, 1>;

#declare KnotsU = OpenUniformKnotVector(Order, NrU, MinUV.u, MaxUV.u)
#declare KnotsV = OpenUniformKnotVector(Order, NrV, MinUV.v, MaxUV.v)

#declare xFn =
  NURBS_SurfFunction(Order, KnotsU, KnotsV, CtrlPoints, Weights, x)
#declare yFn =
  NURBS_SurfFunction(Order, KnotsU, KnotsV, CtrlPoints, Weights, y)
#declare zFn =
  NURBS_SurfFunction(Order, KnotsU, KnotsV, CtrlPoints, Weights, z)

object {
  ParametricMesh(
    function(u, v) { xFn(u, v) },
    function(u, v) { yFn(u, v) },
    function(u, v) { zFn(u, v) },
    MinUV, // Also try <0.4, 0>
    MaxUV, // Also try <0.8, 1>
    <50, 50> // Also try <5, 5>
  )
  pigment { color Blue + Cyan/3 }
}

object {
  CylinderMesh(
    function(u, v) { xFn(u, v) },
    function(u, v) { yFn(u, v) },
    function(u, v) { zFn(u, v) },
    MinUV, // Also try <0.4, 0>
    MaxUV, // Also try <0.8, 1>
    <50, 50> // Also try <5, 5>
    0.02
  )
  pigment { color Magenta/2 + White/4 + Red/2 }
}
*/

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// The photo gear

light_source { <2, 3, -5>*10 color White shadowless }

camera {
  location <1, 1, -2>*10
  look_at <-1, -7, 0>
}

background { color Blue/4 }

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7


Post a reply to this message

Copyright 2003-2023 Persistence of Vision Raytracer Pty. Ltd.