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
<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
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 */
- define N_SEG 20
- define plot(x, y) put_pixel_clip(img, x, y, r, g, b)
- 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; }
- if 0
/* draw only points */ for (i=0; i <= N_SEG; ++i) { plot( pts[i][0], pts[i][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] ); }
- endif
}
- undef plot
- undef line</lang>
Fortran
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>
PHP
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>
Python
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
<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
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
<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>