Bitmap/Bézier curves/Cubic

From Rosetta Code
Revision as of 13:51, 7 November 2009 by rosettacode>Nico80 (Added R code)
Task
Bitmap/Bézier curves/Cubic
You are encouraged to solve this task according to the task description, using any language you may know.

Using the data storage type defined on this page for raster images, and the draw_line function defined in this other one, draw a cubic bezier curves (definition on Wikipedia).


Ada

<lang ada> procedure Cubic_Bezier

         (  Picture        : in out Image;
            P1, P2, P3, P4 : Point;
            Color          : Pixel;
            N              : Positive := 20
         )  is
  Points : array (0..N) of Point;

begin

  for I in Points'Range loop
     declare
        T : constant Float := Float (I) / Float (N);
        A : constant Float := (1.0 - T)**3;
        B : constant Float := 3.0 * T * (1.0 - T)**2;
        C : constant Float := 3.0 * T**2 * (1.0 - T);
        D : constant Float := T**3;
     begin
        Points (I).X := Positive (A * Float (P1.X) + B * Float (P2.X) + C * Float (P3.X) + D * Float (P4.X));
        Points (I).Y := Positive (A * Float (P1.Y) + B * Float (P2.Y) + C * Float (P3.Y) + D * Float (P4.Y));
     end;
  end loop;
  for I in Points'First..Points'Last - 1 loop
     Line (Picture, Points (I), Points (I + 1), Color);
  end loop;

end Cubic_Bezier; </lang> The following test <lang ada>

  X : Image (1..16, 1..16);

begin

  Fill (X, White);
  Cubic_Bezier (X, (16, 1), (1, 4), (3, 16), (15, 11), Black);
  Print (X);

</lang> should produce output:





       HH
     HH  HH
    H      H
    H      H
   H       H
  H        H
 H         H
 H         H
 H         H
 H         H
H         H
H

ALGOL 68

Translation of: Ada
Works with: ALGOL 68 version Standard - pragmat read is an extension
Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386

<lang algol>PRAGMAT READ "Bresenhams_line_algorithm.a68" PRAGMAT;

cubic bezier OF class image :=

         (  REF IMAGE picture,
            POINT p1, p2, p3, p4,
            PIXEL color,
            UNION(INT, VOID) in n
         )VOID:

BEGIN

  INT n = (in n|(INT n):n|20); # default 20 #
  [0:n]POINT points;
  FOR i FROM LWB points TO UPB points DO
        REAL t = i / n,
             a = (1 - t)**3,
             b = 3 * t * (1 - t)**2,
             c = 3 * t**2 * (1 - t),
             d = t**3;
        x OF points [i] := ENTIER (0.5 + a * x OF p1 + b * x OF p2 + c * x OF p3 + d * x OF p4);
        y OF points [i] := ENTIER (0.5 + a * y OF p1 + b * y OF p2 + c * y OF p3 + d * y OF p4)
  OD;
  FOR i FROM LWB points TO UPB points - 1 DO
     (line OF class image)(picture, points (i), points (i + 1), color)
  OD

END # cubic bezier #;

The following test

IF test THEN

  REF IMAGE x = INIT LOC[16,16]PIXEL;
  (fill OF class image)(x, (white OF class image));
  (cubic bezier OF class image)(x, (16, 1), (1, 4), (3, 16), (15, 11), (black OF class image), EMPTY);
  (print OF class image) (x)

FI</lang> Output:

ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
ffffffffffffffffffffffffffffffffffffffffff000000000000ffffffffffffffffffffffffffffffffffffffffff
ffffffffffffffffffffffffffffff000000000000ffffffffffff000000000000ffffffffffffffffffffffffffffff
ffffffffffffffffffffffff000000ffffffffffffffffffffffffffffffffffff000000ffffffffffffffffffffffff
ffffffffffffffffffffffff000000ffffffffffffffffffffffffffffffffffff000000ffffffffffffffffffffffff
ffffffffffffffffff000000ffffffffffffffffffffffffffffffffffffffffff000000ffffffffffffffffffffffff
ffffffffffff000000ffffffffffffffffffffffffffffffffffffffffffffffff000000ffffffffffffffffffffffff
ffffff000000ffffffffffffffffffffffffffffffffffffffffffffffffffffff000000ffffffffffffffffffffffff
ffffff000000ffffffffffffffffffffffffffffffffffffffffffffffffffffff000000ffffffffffffffffffffffff
ffffff000000ffffffffffffffffffffffffffffffffffffffffffffffffffffff000000ffffffffffffffffffffffff
ffffff000000ffffffffffffffffffffffffffffffffffffffffffffffffffffff000000ffffffffffffffffffffffff
000000ffffffffffffffffffffffffffffffffffffffffffffffffffffff000000ffffffffffffffffffffffffffffff
000000ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff

C

"Interface" imglib.h.

<lang c>void cubic_bezier(

      	image img,
       unsigned int x1, unsigned int y1,
       unsigned int x2, unsigned int y2,
       unsigned int x3, unsigned int y3,
       unsigned int x4, unsigned int y4,
       color_component r,
       color_component g,
       color_component b );</lang>

<lang c>#include <math.h>

/* number of segments for the curve */

  1. define N_SEG 20
  1. define plot(x, y) put_pixel_clip(img, x, y, r, g, b)
  2. define line(x0,y0,x1,y1) draw_line(img, x0,y0,x1,y1, r,g,b)

void cubic_bezier(

      	image img,
       unsigned int x1, unsigned int y1,
       unsigned int x2, unsigned int y2,
       unsigned int x3, unsigned int y3,
       unsigned int x4, unsigned int y4,
       color_component r,
       color_component g,
       color_component b )

{

   unsigned int i;
   double pts[N_SEG+1][2];
   for (i=0; i <= N_SEG; ++i)
   {
       double t = (double)i / (double)N_SEG;
       double a = pow((1.0 - t), 3.0);
       double b = 3.0 * t * pow((1.0 - t), 2.0);
       double c = 3.0 * pow(t, 2.0) * (1.0 - t);
       double d = pow(t, 3.0);
       double x = a * x1 + b * x2 + c * x3 + d * x4;
       double y = a * y1 + b * y2 + c * y3 + d * y4;
       pts[i][0] = x;
       pts[i][1] = y;
   }

  1. if 0
   /* draw only points */
   for (i=0; i <= N_SEG; ++i)
   {
       plot( pts[i][0],
             pts[i][1] );
   }
  1. else
   /* draw segments */
   for (i=0; i < N_SEG; ++i)
   {
       int j = i + 1;

line( pts[i][0], pts[i][1],

             pts[j][0], pts[j][1] );
   }
  1. endif

}

  1. undef plot
  2. undef line</lang>

Fortran

Translation of: C

This subroutine should go inside the RCImagePrimitive module (see Bresenham's line algorithm)

<lang fortran> subroutine cubic_bezier(img, p1, p2, p3, p4, color)

   type(rgbimage), intent(inout) :: img
   type(point), intent(in) :: p1, p2, p3, p4
   type(rgb), intent(in) :: color
   integer :: i, j
   real :: pts(0:N_SEG,0:1), t, a, b, c, d, x, y
   do i = 0, N_SEG
      t = real(i) / real(N_SEG)
      a = (1.0 - t)**3.0
      b = 3.0 * t * (1.0 - t)**2.0
      c = 3.0 * (1.0 - t) * t**2.0
      d = t**3.0
      x = a * p1%x + b * p2%x + c * p3%x + d * p4%x
      y = a * p1%y + b * p2%y + c * p3%y + d * p4%y
      pts(i,0) = x
      pts(i,1) = y
   end do
   do i = 0, N_SEG-1
      j = i + 1
      call draw_line(img, point(pts(i,0), pts(i,1)), &
                     point(pts(j,0), pts(j,1)), color)
   end do
 end subroutine cubic_bezier</lang>

J

Solution:
See the Bernstein Polynomials essay on the J Wiki.
Uses code from Basic bitmap storage, Bresenham's line algorithm and Midpoint circle algorithm. <lang j>require 'numeric'

bik=: 2 : '((*&(u!v))@(^&u * ^&(v-u)@-.))' basiscoeffs=: <: 4 : 'x bik y t. i.>:y'"0~ i. linearcomb=: basiscoeffs@#@[ evalBernstein=: ([ +/ .* linearcomb) p. ] NB. evaluate Bernstein Polynomial (general)

NB.*getBezierPoints v Returns points for bezier curve given control points (y) NB. eg: getBezierPoints controlpoints NB. y is: y0 x0, y1 x1, y2 x2 ... getBezierPoints=: monad define

 ctrlpts=. (/: {:"1)  _2]\ y  NB. sort ctrlpts for increasing x
 xvals=. ({: ,~ {. + +:@:i.@<.@-:@-~/) ({:"1) 0 _1{ctrlpts
 tvals=.  ((] - {.) % ({: - {.)) xvals
 xvals ,.~ ({."1 ctrlpts) evalBernstein tvals

)

NB.*drawBezier v Draws bezier curve defined by (x) on image (y) NB. eg: (42 40 10 30 186 269 26 187;255 0 0) drawBezier myimg NB. x is: 2-item list of boxed (controlpoints) ; (color) drawBezier=: (1&{:: ;~ 2 ]\ [: roundint@getBezierPoints"1 (0&{::))@[ drawLines ]</lang>

Example usage: <lang j> myimg=: 0 0 255 makeRGB 300 300

  ]randomctrlpts=: ,3 2 ?@$ }:$ myimg                               NB. 3 control points - quadratic
  ]randomctrlpts=: ,4 2 ?@$ }:$ myimg                               NB. 4 control points - cubic
  myimg=: ((2 ,.~ _2]\randomctrlpts);255 0 255) drawCircles myimg   NB. draw control points
  viewRGB (randomctrlpts; 255 255 0) drawBezier myimg               NB. display image with bezier line</lang>

OCaml

<lang ocaml>let cubic_bezier ~img ~color

       ~p1:(_x1, _y1)
       ~p2:(_x2, _y2)
       ~p3:(_x3, _y3)
       ~p4:(_x4, _y4) =
 let x1, y1, x2, y2, x3, y3, x4, y4 =
   (float _x1, float _y1,
    float _x2, float _y2,
    float _x3, float _y3,
    float _x4, float _y4)
 in
 let bz t =
   let a = (1.0 -. t) ** 3.0
   and b = 3.0 *. t *. ((1.0 -. t) ** 2.0)
   and c = 3.0 *. (t ** 2.0) *. (1.0 -. t)
   and d = t ** 3.0
   in
   let x = a *. x1 +. b *. x2 +. c *. x3 +. d *. x4
   and y = a *. y1 +. b *. y2 +. c *. y3 +. d *. y4
   in
   (int_of_float x, int_of_float y)
 in
 let rec loop _t acc =
   if _t > 20 then acc else
   begin
     let t = (float _t) /. 20.0 in
     let x, y = bz t in
     loop (succ _t) ((x,y)::acc)
   end
 in
 let pts = loop 0 [] in
 (*
 (* draw only points *)
 List.iter (fun (x, y) -> put_pixel img color x y) pts;
 *)
 (* draw segments *)
 let line = draw_line ~img ~color in
 let by_pair li f =
   let rec aux prev = function
     | [] -> ()
     | x::xs ->
         f prev x;
         aux x xs
   in
   aux (List.hd li) (List.tl li)
 in
 by_pair pts (fun p0 p1 -> line ~p0 ~p1);
</lang>

R

<lang R>

  1. x, y: the x and y coordinates of the hull points
  2. n: the number of points in the curve.

bezierCurve <- function(x, y, n=10) { outx <- NULL outy <- NULL

i <- 1 for (t in seq(0, 1, length.out=n)) { b <- bez(hullx, hully, t) outx[i] <- b$x outy[i] <- b$y

i <- i+1 }

return (list(x=outx, y=outy)) }

bez <- function(x, y, t) { outx <- 0 outy <- 0 n <- length(hullx)-1 for (i in 0:n) { outx <- outx + choose(n, i)*((1-t)^(n-i))*t^i*hullx[i+1] outy <- outy + choose(n, i)*((1-t)^(n-i))*t^i*hully[i+1] }

return (list(x=outx, y=outy)) }

  1. Example usage

x <- c(4,6,4,5,6,7) y <- 1:6 plot(x, y, "o", pch=20) points(bezierCurve(x,y,20), type="l", col="red") </lang>

Ruby

Translation of: Tcl

<lang ruby>class Pixmap

 def draw_bezier_curve(points, colour)
   # ensure the points are increasing along the x-axis
   points = points.sort_by {|p| [p.x, p.y]}
   xmin = points[0].x
   xmax = points[-1].x
   increment = 2
   prev = points[0]
   ((xmin + increment) .. xmax).step(increment) do |x|
     t = 1.0 * (x - xmin) / (xmax - xmin)
     p = Pixel[x, bezier(t, points).round]
     draw_line(prev, p, colour)
     prev = p
   end
 end

end

  1. the generalized n-degree Bezier summation

def bezier(t, points)

 n = points.length - 1
 points.each_with_index.inject(0.0) do |sum, (point, i)|
   sum += n.choose(i) * (1-t)**(n - i) * t**i * point.y
 end

end

class Fixnum

 def choose(k)
   self.factorial / (k.factorial * (self - k).factorial)
 end
 def factorial
   (2 .. self).reduce(1, :*)
 end

end

bitmap = Pixmap.new(400, 400) points = [

 Pixel[40,100], Pixel[100,350], Pixel[150,50], 
 Pixel[150,150], Pixel[350,250], Pixel[250,250]

] points.each {|p| bitmap.draw_circle(p, 3, RGBColour::RED)} bitmap.draw_bezier_curve(points, RGBColour::BLUE)</lang>

Tcl

Library: Tk

This solution can be applied to any number of points. Uses code from Basic bitmap storage (newImage, fill), Bresenham's line algorithm (drawLine), and Midpoint circle algorithm (drawCircle) <lang tcl>package require Tcl 8.5 package require Tk

proc drawBezier {img colour args} {

   # ensure the points are increasing along the x-axis
   set points [lsort -real -index 0 $args]
   
   set xmin [x [lindex $points 0]]
   set xmax [x [lindex $points end]]
   set prev [lindex $points 0]
   set increment 2
   for {set x [expr {$xmin + $increment}]} {$x <= $xmax} {incr x $increment} {
       set t [expr {1.0 * ($x - $xmin) / ($xmax - $xmin)}]
       set this [list $x [::tcl::mathfunc::round [bezier $t $points]]]
       drawLine $img $colour $prev $this
       set prev $this
   }

}

  1. the generalized n-degree Bezier summation

proc bezier {t points} {

   set n [expr {[llength $points] - 1}]
   for {set i 0; set sum 0.0} {$i <= $n} {incr i} {
       set sum [expr {$sum + [C $n $i] * (1-$t)**($n - $i) * $t**$i * [y [lindex $points $i]]}]
   }
   return $sum

}

proc C {n i} {expr {[ifact $n] / ([ifact $i] * [ifact [expr {$n - $i}]])}} proc ifact n {

   for {set i $n; set sum 1} {$i >= 2} {incr i -1} {
       set sum [expr {$sum * $i}]
   }
   return $sum

}

proc x p {lindex $p 0} proc y p {lindex $p 1}

proc newbezier {n w} {

   set size 400
   set bezier [newImage $size $size]
   fill $bezier white
   for {set i 1} {$i <= $n} {incr i} {
       set point [list [expr {int($size*rand())}] [expr {int($size*rand())}]]
       lappend points $point
       drawCircle $bezier red $point 3
   }
   puts $points
   
   drawBezier $bezier blue {*}$points    
   
   $w configure -image $bezier

}

set degree 4 ;# cubic bezier -- for quadratic, use 3 label .img button .new -command [list newbezier $degree .img] -text New button .exit -command exit -text Exit pack .new .img .exit -side top</lang> Results in:

TI-89 BASIC

Note: This example does not use a user-defined image type, since that would be particularly impractical, but rather draws on the calculator's graph screen, which has essentially the same operations as an implementation of Basic bitmap storage would, except for being black-and-white.
Define cubic(p1,p2,p3,p4,segs) = Prgm
  Local i,t,u,prev,pt
  0 → pt
  For i,1,segs+1
    (i-1.0)/segs → t   © Decimal to avoid slow exact arithetic
    (1-t) → u
    pt → prev
    u^3*p1 + 3t*u^2*p2 + 3t^2*u*p3 + t^3*p4 → pt
    If i>1 Then
      PxlLine floor(prev[1,1]), floor(prev[1,2]), floor(pt[1,1]), floor(pt[1,2])
    EndIf
  EndFor
EndPrgm