Death Star

From Rosetta Code
Revision as of 07:49, 21 August 2011 by 79.54.59.51 (talk) (Added Python version)
Task
Death Star
You are encouraged to solve this task according to the task description, using any language you may know.

Death Star is a task to display a region that consists of a large sphere with part of a smaller sphere removed from it as a result of geometric subtraction. (This will basically produce a shape like a "death star".)

See also: Draw a sphere.

Brlcad

<lang brlcad># We need a database to hold the objects opendb deathstar.g y

  1. We will be measuring in kilometers

units km

  1. Create a sphere of radius 60km centred at the origin

in sph1.s sph 0 0 0 60

  1. We will be subtracting an overlapping sphere with a radius of 40km
  2. The resultant hole will be smaller than this, because we only
  3. only catch the edge

in sph2.s sph 0 90 0 40

  1. Create a region named deathstar.r which consists of big minus small sphere

r deathstar.r u sph1.s - sph2.s

  1. We will use a plastic material texture with rgb colour 224,224,224
  2. with specular lighting value of 0.1 and no inheritance

mater deathstar.r "plastic sp=0.1" 224 224 224 0

  1. Clear the wireframe display and draw the deathstar

B deathstar.r

  1. We now trigger the raytracer to see our finished product

rt</lang>

C

Primitive ray tracing. <lang c>#include <stdio.h>

  1. include <math.h>

char shades[] = ".:!*oe&#%@";

double light[3] = { -50, 30, 50 }; void normalize(double * v) { double len = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); v[0] /= len; v[1] /= len; v[2] /= len; }

double dot(double *x, double *y) { double d = x[0]*y[0] + x[1]*y[1] + x[2]*y[2]; return d < 0 ? -d : 0; }

typedef struct { double cx, cy, cz, r; } sphere_t;

/* positive shpere and negative sphere */ sphere_t pos = { 20, 20, 0, 20 }, neg = { 1, 1, -6, 20 };

/* check if a ray (x,y, -inf)->(x, y, inf) hits a sphere; if so, return

  the intersecting z values.  z1 is closer to the eye */

int hit_sphere(sphere_t *sph, double x, double y, double *z1, double *z2) { double zsq; x -= sph->cx; y -= sph->cy; zsq = sph->r * sph->r - (x * x + y * y); if (zsq < 0) return 0; zsq = sqrt(zsq); *z1 = sph->cz - zsq; *z2 = sph->cz + zsq; return 1; }

void draw_sphere(double k, double ambient) { int i, j, intensity, hit_result; double b; double vec[3], x, y, zb1, zb2, zs1, zs2; for (i = floor(pos.cy - pos.r); i <= ceil(pos.cy + pos.r); i++) { y = i + .5; for (j = floor(pos.cx - 2 * pos.r); j <= ceil(pos.cx + 2 * pos.r); j++) { x = (j - pos.cx) / 2. + .5 + pos.cx;

/* ray lands in blank space, draw bg */ if (!hit_sphere(&pos, x, y, &zb1, &zb2)) hit_result = 0;

/* ray hits pos sphere but not neg, draw pos sphere surface */ else if (!hit_sphere(&neg, x, y, &zs1, &zs2)) hit_result = 1;

/* ray hits both, but pos front surface is closer */ else if (zs1 > zb1) hit_result = 1;

/* pos sphere surface is inside neg sphere, show bg */ else if (zs2 > zb2) hit_result = 0;

/* back surface on neg sphere is inside pos sphere, the only place where neg sphere surface will be shown */ else if (zs2 > zb1) hit_result = 2; else hit_result = 1;

switch(hit_result) { case 0: putchar(' '); continue; case 1: vec[0] = x - pos.cx; vec[1] = y - pos.cy; vec[2] = zb1 - pos.cz; break; default: vec[0] = neg.cx - x; vec[1] = neg.cy - y; vec[2] = neg.cz - zs2; }

normalize(vec); b = pow(dot(light, vec), k) + ambient; intensity = (1 - b) * (sizeof(shades) - 1); if (intensity < 0) intensity = 0; if (intensity >= sizeof(shades) - 1) intensity = sizeof(shades) - 2; putchar(shades[intensity]); } putchar('\n'); } }

int main() { normalize(light);

draw_sphere(2, .5); return 0; }</lang>output<lang> eeeee:::::::

                               eeeeeeeee..............                          
                            ooeeeeeeeeee..................                      
                          ooooeeeeeeeee......................                   
                       oooooooeeeeeeee..........................                
                     ooooooooooeeeee..............................              
                   **ooooooooooeeee.................................            
                 ****ooooooooooee.....................................          
               !*****ooooooooooe.......................................         
             !!!*****ooooooooo:..........................................       
           :!!!!*****ooooooo:::...........................................      
         :::!!!!*****ooooo!:::::...........................................     
       ::::!!!!!*****ooo!!!!::::............................................    
      .::::!!!!*****oo*!!!!!::::............................................    
    ...::::!!!!*********!!!!:::::............................................   
   ...::::!!!!****o*****!!!!!::::............................................   
 ....::::!!!!***ooo******!!!!!::::............................................  
....::::!!!!*ooooooo*****!!!!!:::::...........................................  

...::::!!!!!oooooooooo*****!!!!!:::::..........................................

!!!!eeooooooooooo******!!!!!:::::.........................................

!!!!!eeeeeeeooooooooooo******!!!!!:::::........................................ eeeeeeeeeeeeoooooooooooo******!!!!!:::::....................................... eeeeeeeeeeeeeoooooooooooo******!!!!!!:::::..................................... eeeeeeeeeeeeeeoooooooooooo******!!!!!!:::::....................................

eeeeeeeeeeeeeeoooooooooooo*******!!!!!!:::::.................................   
eeeeeeeeeeeeeeeooooooooooooo******!!!!!!::::::..............................:   
 eeeeeeeeeeeeeeeooooooooooooo*******!!!!!!:::::::..........................:    
 eeeeeeeeeeeeeeeeoooooooooooooo*******!!!!!!!:::::::.....................::!    
  eeeeeeeeeeeeeeeeeooooooooooooo********!!!!!!!:::::::::..............::::!     
   eeeeeeeeeeeeeeeeeoooooooooooooo********!!!!!!!!::::::::::::::::::::::!*      
    eeeeeeeeeeeeeeeeeeooooooooooooooo********!!!!!!!!!!:::::::::::::!!!!*       
      eeeeeeeeeeeeeeeeeoooooooooooooooo**********!!!!!!!!!!!!!!!!!!!!!*         
       eeeeeeeeeeeeeeeeeeooooooooooooooooo************!!!!!!!!!!!!****          
         eeeeeeeeeeeeeeeeeeoooooooooooooooooo**********************o            
           eeeeeeeeeeeeeeeeeeeooooooooooooooooooooo************ooo              
             eeeeeeeeeeeeeeeeeeeeooooooooooooooooooooooooooooooo                
                eeeeeeeeeeeeeeeeeeeeooooooooooooooooooooooooo                   
                   eeeeeeeeeeeeeeeeeeeeeoooooooooooooooooo                      
                       eeeeeeeeeeeeeeeeeeeeeeeeeeeeeee                          
                              eeeeeeeeeeeeeeeee</lang>

D

Translation of: C

<lang d>import std.stdio, std.math, std.numeric, std.algorithm;

/*const*/ struct V3 {

   double[3] v;
   @property V3 normalize() pure nothrow const {
       //immutable double len = sqrt(dotProduct(v, v));
       immutable double len= sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
       // return V3(v[] / len);
       return V3([v[0] / len, v[1] / len, v[2] / len]);
   }
   double dot(const ref V3 y) pure nothrow const {
       immutable double d = dotProduct(v, y.v);
       return d < 0 ? -d : 0;
   }

}


const struct Sphere { double cx, cy, cz, r; }

void drawSphere(in double k, in double ambient, in V3 light)nothrow{

   /** Check if a ray (x,y, -inf).(x, y, inf) hits a sphere; if so,
   return the intersecting z values.  z1 is closer to the eye */
   static bool hitSphere(const ref Sphere sph,
                         in double x0, in double y0,
                         out double z1,
                         out double z2) pure nothrow {
       immutable double x = x0 - sph.cx;
       immutable double y = y0 - sph.cy;
       immutable double zsq = sph.r ^^ 2 - (x ^^ 2 + y ^^ 2);
       if (zsq < 0)
           return false;
       immutable double szsq = sqrt(zsq);
       z1 = sph.cz - szsq;
       z2 = sph.cz + szsq;
       return true;
   }
   enum string shades = ".:!*oe&#%@";
   // positive and negative spheres
   enum pos = Sphere(20, 20, 0, 20);
   enum neg = Sphere(1, 1, -6, 20);
   foreach (int i; cast(int)floor(pos.cy - pos.r) ..
                   cast(int)ceil(pos.cy + pos.r) + 1) {
       immutable double y = i + 0.5;
       jloop:
       foreach (int j; cast(int)floor(pos.cx - 2 * pos.r) ..
                       cast(int)ceil(pos.cx + 2 * pos.r) + 1) {
           immutable double x = (j - pos.cx) / 2.0 + 0.5 + pos.cx;
           enum Hit { background, posSphere, negSphere }
           Hit hitResult;
           double zb1, zs2;
           {
               double zb2, zs1;
               if (!hitSphere(pos, x, y, zb1, zb2)) {
                   // Ray lands in blank space, draw bg
                   hitResult = Hit.background;
               } else if (!hitSphere(neg, x, y, zs1, zs2)) {
                   // Ray hits pos sphere but not neg one,
                   // draw pos sphere surface
                   hitResult = Hit.posSphere;
               } else if (zs1 > zb1) {
                   // ray hits both, but pos front surface is closer
                   hitResult = Hit.posSphere;
               } else if (zs2 > zb2) {
                   // pos sphere surface is inside neg sphere,
                   // show bg
                   hitResult = Hit.background;
               } else if (zs2 > zb1) {
                   // Back surface on neg sphere is inside pos
                   // sphere, the only place where neg sphere
                   // surface will be shown
                   hitResult = Hit.negSphere;
               } else {
                   hitResult = Hit.posSphere;
               }
           }
           V3 vec_;
           final switch (hitResult) {
               case Hit.background:
                   putchar(' ');
                   continue jloop;
               case Hit.posSphere:
                   vec_ = V3([x - pos.cx, y - pos.cy, zb1 - pos.cz]);
                   break;
               case Hit.negSphere:
                   vec_ = V3([neg.cx-x, neg.cy-y, neg.cz-zs2]);
                   break;
           }
           immutable V3 nvec = vec_.normalize;
           immutable double b = light.dot(nvec) ^^ k + ambient;
           int intensity = cast(int)((1 - b) * shades.length);
           intensity = min(shades.length, max(0, intensity));
           putchar(shades[intensity]);
       }
       putchar('\n');
   }

}


void main() {

   enum light = V3([-50, 30, 50]).normalize;
   drawSphere(2, 0.5, light);

}</lang> The output is the same of the C version.

DWScript

Translation of: C

<lang delphi>const cShades = '.:!*oe&#%@';

type TVector = array [0..2] of Float;

var light : TVector = [-50.0, 30, 50];

procedure Normalize(var v : TVector); begin

  var len := Sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
  v[0] /= len; v[1] /= len; v[2] /= len;

end;

function Dot(x, y : TVector) : Float; begin

  var d :=x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
  if d<0 then 
     Result:=-d 
  else Result:=0;

end;

type

  TSphere = record
     cx, cy, cz, r : Float;
  end;
  

const big : TSphere = (cx: 20; cy: 20; cz: 0; r: 20); const small : TSphere = (cx: 7; cy: 7; cz: -10; r: 15);

function HitSphere(sph : TSphere; x, y : Float; var z1, z2 : Float) : Boolean; begin

  x -= sph.cx;
  y -= sph.cy;
  var zsq = sph.r * sph.r - (x * x + y * y);
  if (zsq < 0) then Exit False;
  zsq := Sqrt(zsq);
  z1 := sph.cz - zsq;
  z2 := sph.cz + zsq;
  Result:=True;

end;

procedure DrawSphere(k, ambient : Float); var

  i, j, intensity : Integer;
  b : Float;
  x, y, zb1, zb2, zs1, zs2 : Float;
  vec : TVector;

begin

  for i:=Trunc(big.cy-big.r) to Trunc(big.cy+big.r)+1 do begin
     y := i + 0.5;
     for j := Trunc(big.cx-2*big.r) to Trunc(big.cx+2*big.r) do begin
        x := (j-big.cx)/2 + 0.5 + big.cx;

        if not HitSphere(big, x, y, zb1, zb2) then begin
           Print(' ');
           continue;
        end;
        if not HitSphere(small, x, y, zs1, zs2) then begin
           vec[0] := x - big.cx;
           vec[1] := y - big.cy;
           vec[2] := zb1 - big.cz;
        end else begin
           if zs1 < zb1 then begin
              if zs2 > zb2 then begin
                 Print(' ');
                 continue;
              end;
              if zs2 > zb1 then begin
                 vec[0] := small.cx - x;
                 vec[1] := small.cy - y;
                 vec[2] := small.cz - zs2;
              end else begin
                 vec[0] := x - big.cx;
                 vec[1] := y - big.cy;
                 vec[2] := zb1 - big.cz;
              end;
           end else begin
              vec[0] := x - big.cx;
              vec[1] := y - big.cy;
              vec[2] := zb1 - big.cz;
           end;
        end;

        Normalize(vec);
        b := Power(Dot(light, vec), k) + ambient;
        intensity := Round((1 - b) * Length(cShades));
        Print(cShades[ClampInt(intensity+1, 1, Length(cShades))]);
     end;
     PrintLn();
  end;

end;

Normalize(light);

DrawSphere(2, 0.3);</lang>

Perl

Writes a PGM to stdout. <lang perl>use strict;

sub sq { my $s = 0; $s += $_ ** 2 for @_; $s; }

sub hit { my ($sph, $x, $y) = @_; $x -= $sph->[0]; $y -= $sph->[1];

my $z = sq($sph->[3]) - sq($x, $y); return if $z < 0;

$z = sqrt $z; return $sph->[2] - $z, $sph->[2] + $z; }

sub normalize { my $v = shift; my $n = sqrt sq(@$v); $_ /= $n for @$v; $v; }

sub dot { my ($x, $y) = @_; my $s = $x->[0] * $y->[0] + $x->[1] * $y->[1] + $x->[2] * $y->[2]; $s > 0 ? $s : 0; }

my $pos = [ 120, 120, 0, 120 ]; my $neg = [ -77, -33, -100, 190 ]; my $light = normalize([ -12, 13, -10 ]); sub draw { my ($k, $amb) = @_; binmode STDOUT, ":raw"; print "P5\n", $pos->[0] * 2 + 3, " ", $pos->[1] * 2 + 3, "\n255\n"; for my $y (($pos->[1] - $pos->[3] - 1) .. ($pos->[1] + $pos->[3] + 1)) { my @row = (); for my $x (($pos->[0] - $pos->[3] - 1) .. ($pos->[0] + $pos->[3] + 1)) { my ($hit, @hs) = 0; my @h = hit($pos, $x, $y);

if (!@h) { $hit = 0 } elsif (!(@hs = hit($neg, $x, $y))) { $hit = 1 } elsif ($hs[0] > $h[0]) { $hit = 1 } elsif ($hs[1] > $h[0]) { $hit = $hs[1] > $h[1] ? 0 : 2 } else { $hit = 1 }

my ($val, $v); if ($hit == 0) { $val = 0 } elsif ($hit == 1) { $v = [ $x - $pos->[0], $y - $pos->[1], $h[0] - $pos->[2] ]; } else { $v = [ $neg->[0] - $x, $neg->[1] - $y, $neg->[2] - $hs[1] ]; } if ($v) { normalize($v); $val = int((dot($v, $light) ** $k + $amb) * 255); $val = ($val > 255) ? 255 : ($val < 0) ? 0 : $val; } push @row, $val; } print pack("C*", @row); } }

draw(2, 0.2);</lang>

Perl 6

Translation of: C

Reimplemented to output a .pgm image.

<lang perl6>class sphere {

  has $.cx; # center x coordinate
  has $.cy; # center y coordinate
  has $.cz; # center z coordinate
  has $.r;  # radius

}

my $depth = 255; # image color depth

my $x = my $y = 255; # dimensions of generated .pgm; must be odd

my $s = ($x - 1)/2; # scaled dimension to build geometry

my @light = normalize([ -2, 4, -5 ]);

  1. positive sphere at origin

my $pos = sphere.new(

   cx => 0,
   cy => 0,
   cz => 0,
   r  => $s.Int

);

  1. negative sphere offset to upper left

my $neg = sphere.new(

   cx => (-$s*.95).Int,
   cy => (-$s*.95).Int,
   cz => (-$s*.2).Int,
   r  => ($s*.7).Int

);

sub MAIN ($outfile = 'deathstar-perl6.pgm') {

   my $out = open( $outfile, :w, :bin ) or die "$!\n";
   $out.say("P5\n$x $y\n$depth"); # .pgm header
   say 'Calculating row:';
   $out.print( draw_ds(3, .15)».chrs );
   $out.close;

}

sub draw_ds ( $k, $ambient ) {

   my @pixels;
   my $bs = "\b" x 8;
   for ($pos.cy - $pos.r) .. ($pos.cy + $pos.r) -> $y {
       note $bs, $y, ' '; # monitor progress
       for ($pos.cx - $pos.r) .. ($pos.cx + $pos.r) -> $x {
           # black if we don't hit positive sphere, ignore negative sphere
           if not hit($pos, $x, $y, my $posz) {
               @pixels.push(0);
               next;
           }
           my @vec;
           # is front of positive sphere inside negative sphere?
           if hit($neg, $x, $y, my $negz) and $negz.min < $posz.min < $negz.max {
               # make black if whole positive sphere eaten here
               if $negz.min < $posz.max < $negz.max { @pixels.push(0); next; }
               # render inside of negative sphere
               @vec = normalize([$neg.cx - $x, $neg.cy - $y, -$negz.max - $neg.cz]);
           }
           else {
               # render outside of positive sphere
               @vec = normalize([$x - $pos.cx, $y - $pos.cy,  $posz.max - $pos.cz]);
           }
           my $intensity = dot(@light, @vec) ** $k + $ambient;
           @pixels.push( ($intensity * $depth).Int min $depth );
       }
   }
   say $bs, 'Writing file.';
   return @pixels;

}

  1. normalize a vector

sub normalize (@vec) { return @vec »/» ([+] @vec Z* @vec).sqrt }

  1. dot product of two vectors

sub dot (@x, @y) { return -([+] @x Z* @y) max 0 }

  1. are the coordinates within the radius of the sphere?

sub hit ($sphere, $x is copy, $y is copy, $z is rw) {

   $x -= $sphere.cx;
   $y -= $sphere.cy;
   my $z2 = $sphere.r * $sphere.r - $x * $x - $y * $y;
   return 0 if $z2 < 0;
   $z2 = $z2.sqrt;
   $z = $sphere.cz - $z2 .. $sphere.cz + $z2;
   return 1;

}</lang>

Python

Translation of: C

<lang python>import sys, math, collections

Sphere = collections.namedtuple("Sphere", "cx cy cz r") V3 = collections.namedtuple("V3", "x y z")

def normalize((x, y, z)):

   len = math.sqrt(x**2 + y**2 + z**2)
   return V3(x / len, y / len, z / len)

def dot(v1, v2):

   d = v1.x*v2.x + v1.y*v2.y + v1.z*v2.z
   return -d if d < 0 else 0.0

def hit_sphere(sph, x0, y0):

   x = x0 - sph.cx
   y = y0 - sph.cy
   zsq = sph.r ** 2 - (x ** 2 + y ** 2)
   if zsq < 0:
       return (False, 0, 0)
   szsq = math.sqrt(zsq)
   return (True, sph.cz - szsq, sph.cz + szsq)

def draw_sphere(k, ambient, light):

   shades = ".:!*oe&#%@"
   pos = Sphere(20.0, 20.0, 0.0, 20.0)
   neg = Sphere(1.0, 1.0, -6.0, 20.0)
   for i in xrange(int(math.floor(pos.cy - pos.r)),
                   int(math.ceil(pos.cy + pos.r) + 1)):
       y = i + 0.5
       for j in xrange(int(math.floor(pos.cx - 2 * pos.r)),
                       int(math.ceil(pos.cx + 2 * pos.r) + 1)):
           x = (j - pos.cx) / 2.0 + 0.5 + pos.cx
           (h, zb1, zb2) = hit_sphere(pos, x, y)
           if not h:
               hit_result = 0
           else:
               (h, zs1, zs2) = hit_sphere(neg, x, y)
               if not h:
                   hit_result = 1
               elif zs1 > zb1:
                   hit_result = 1
               elif zs2 > zb2:
                   hit_result = 0
               elif zs2 > zb1:
                   hit_result = 2
               else:
                   hit_result = 1
           if hit_result == 0:
               sys.stdout.write(' ')
               continue
           elif hit_result == 1:
               vec = V3(x - pos.cx, y - pos.cy, zb1 - pos.cz)
           elif hit_result == 2:
               vec = V3(neg.cx-x, neg.cy-y, neg.cz-zs2)
           vec = normalize(vec)
           b = dot(light, vec) ** k + ambient
           intensity = int((1 - b) * len(shades))
           intensity = min(len(shades), max(0, intensity))
           sys.stdout.write(shades[intensity])
       print

light = normalize(V3(-50, 30, 50)) draw_sphere(2, 0.5, light)</lang>

Openscad

<lang openscad>// We are performing geometric subtraction

difference() {

 // Create the primary sphere of radius 60 centred at the origin
 translate(v = [0,0,0]) {
   sphere(60);
 }
 /*Subtract an overlapping sphere with a radius of 40
    The resultant hole will be smaller than this, because we only
    only catch the edge
 */
 translate(v = [0,90,0]) {
   sphere(40);
 }

}</lang>

POV-Ray

<lang POV-Ray>camera { perspective location <0.0 , .8 ,-3.0> look_at 0

        aperture .1 blur_samples 20 variance 1/100000 focal_point 0}
                           

light_source{< 3,3,-3> color rgb 1}

sky_sphere { pigment{ color rgb <0,.2,.5>}}

plane {y,-5 pigment {color rgb .54} normal {hexagon} }

difference {

sphere { 0,1 }
sphere { <-1,1,-1>,1 }
 texture { 
   pigment{ granite } 
   finish { phong 1 reflection {0.10 metallic 0.5} }
 } 

} </lang>