Bitmap/Bézier curves/Cubic

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 algol68>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

Factor

The points should probably be in a sequence... <lang factor>USING: arrays kernel locals math math.functions

rosettacode.raster.storage sequences ;

IN: rosettacode.raster.line

! this gives a function

(cubic-bezier) ( P0 P1 P2 P3 -- bezier )

   [ :> x
       1 x - 3 ^ P0 n*v
       1 x - sq 3 * x * P1 n*v
       1 x - 3 * x sq * P2 n*v
       x 3 ^ P3 n*v
       v+ v+ v+ ] ; inline

! gives an interval of x from 0 to 1 to map the bezier function

t-interval ( x -- interval )

   [ iota ] keep 1 - [ / ] curry map ;

! turns a list of points into the list of lines between them

points-to-lines ( seq -- seq )

   dup rest [ 2array ] 2map ;

draw-lines ( {R,G,B} points image -- )

   [ [ first2 ] dip draw-line ] curry with each ;

bezier-lines ( {R,G,B} P0 P1 P2 P3 image -- )

   ! 100 is an arbitrary value.. could be given as a parameter..
   100 t-interval P0 P1 P2 P3 (cubic-bezier) map
   points-to-lines
   {R,G,B} swap image draw-lines ;</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>

PHP

Translation of: Python

Works with: PHP version 4.3.0

Library: GD

Outputs image to the right directly to browser or stdout.

<lang php><?

$image = imagecreate(200, 200); // The first allocated color will be the background color: imagecolorallocate($image, 255, 255, 255); $color = imagecolorallocate($image, 255, 0, 0); cubicbezier($image, $color, 160, 10, 10, 40, 30, 160, 150, 110); imagepng($image);

function cubicbezier($img, $col, $x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3, $n = 20) { $pts = array();

for($i = 0; $i <= $n; $i++) { $t = $i / $n; $t1 = 1 - $t; $a = pow($t1, 3); $b = 3 * $t * pow($t1, 2); $c = 3 * pow($t, 2) * $t1; $d = pow($t, 3);

$x = round($a * $x0 + $b * $x1 + $c * $x2 + $d * $x3); $y = round($a * $y0 + $b * $y1 + $c * $y2 + $d * $y3); $pts[$i] = array($x, $y); }

for($i = 0; $i < $n; $i++) { imageline($img, $pts[$i][0], $pts[$i][1], $pts[$i+1][0], $pts[$i+1][1], $col); } } </lang>

PureBasic

<lang PureBasic>Procedure cubic_bezier(img, p1x, p1y, p2x, p2y, p3x, p3y, p4x, p4y, Color, n_seg)

 Protected i
 Protected.f t, t1, a, b, c, d
 Dim pts.POINT(n_seg)
 
 For i = 0 To n_seg
   t = i / n_seg
   t1 = 1.0 - t
   a = Pow(t1, 3)
   b = 3.0 * t * Pow(t1, 2)
   c = 3.0 * Pow(t, 2) * t1
   d = Pow(t, 3)
   pts(i)\x = a * p1x + b * p2x + c * p3x + d * p4x
   pts(i)\y = a * p1y + b * p2y + c * p3y + d * p4y
 Next
 
 StartDrawing(ImageOutput(img))
   FrontColor(Color)
   For i = 0 To n_seg - 1
     BresenhamLine(pts(i)\x, pts(i)\y, pts(i + 1)\x, pts(i + 1)\y) ;this calls the implementation of a draw_line routine
   Next
 StopDrawing()

EndProcedure

Define w, h, img w = 200: h = 200: img = 1 CreateImage(img, w, h) ;img is internal id of the image

OpenWindow(0, 0, 0, w, h,"Bezier curve, cubic", #PB_Window_SystemMenu) cubic_bezier(1, 160,10, 10,40, 30,160, 150,110, RGB(255, 255, 255), 20) ImageGadget(0, 0, 0, w, h, ImageID(1))

Define event Repeat

 event = WaitWindowEvent()

Until event = #PB_Event_CloseWindow</lang>

Python

Works with: Python version 3.1

Extending the example given here and using the algorithm from the C solution above: <lang python>def cubicbezier(self, x0, y0, x1, y1, x2, y2, x3, y3, n=20):

   pts = []
   for i in range(n+1):
       t = i / n
       a = (1. - t)**3
       b = 3. * t * (1. - t)**2
       c = 3.0 * t**2 * (1.0 - t)
       d = t**3
       
       x = int(a * x0 + b * x1 + c * x2 + d * x3)
       y = int(a * y0 + b * y1 + c * y2 + d * y3)
       pts.add( (x, y) )
   for i in range(n):
       self.line(pts[i][0], pts[i][1], pts[i+1][0], pts[i+1][1])

Bitmap.cubicbezier = cubicbezier

bitmap = Bitmap(17,17) bitmap.cubicbezier(16,1, 1,4, 3,16, 15,11) bitmap.chardisplay()

The origin, 0,0; is the lower left, with x increasing to the right, and Y increasing upwards.

The chardisplay above produces the following output : +-----------------+ | | | | | | | | | @@@@ | | @@@ @@@ | | @ | | @ | | @ | | @ | | @ | | @ | | @ | | @ | | @@@@ | | @@@@| | | +-----------------+ </lang>

R

<lang R># x, y: the x and y coordinates of the hull points

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(x, y, 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(x)-1 for (i in 0:n) { outx <- outx + choose(n, i)*((1-t)^(n-i))*t^i*x[i+1] outy <- outy + choose(n, i)*((1-t)^(n-i))*t^i*y[i+1] }

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

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

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
   }

}

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.

<lang ti89b>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</lang>