// Last updated 2008/12/14 09:00

/*
Originally inspired by:
http://www.imagemagick.org/discourse-server/viewtopic.php?f=2&t=12530
The idea for these specific examples came from reading this:
http://www.csl.mtu.edu/cs4611/www/HillLectureNotes/CS4611%202D%20Affine%20Transformation.htm
When reading that (and other web pages about affine) keep in mind
that IM's ordering of the affine matrix as described at:
http://imagemagick.org/script/command-line-options.php#affine
orders the affine values and their multiplication like this:
[x y 1] |sx rx 0|
|ry sy 0|
|tx ty 1|

Whereas the CS4611 web page uses this (which, if nothing else, is tidier):
|sx ry tx| |x|
|rx sy ty| |y|
|0  0  1 | |1|

My multiplication routine is written to conform to the way IM
specifies things.

ALSO, I think there are a couple of errors on the CS4611 page.
1. In the example of rotation about a point, it says that first
translate by V, then rotate, then translate by -V.
But the matrix representation of this does -V,rotate,V.
2. Reflection across the x-axis is not correct as shown.
When a point (x,y) is reflected across the x-axis its
new coordinate is (x,-y) - the matrix shown in the example
actually reflects across the y-axis - i.e. it produces (-x,y).
*/
#include <windows.h>
#include <wand/magick_wand.h>
#define _USE_MATH_DEFINES
#include <math.h>
#define DegreesToRadians(a) (a*M_PI/180.)
#define RadiansToDegrees(a) (a*180./M_PI)
// Each of these affine functions assumes that the input arrays
// have been correctly defined to be arrays of 6 (or more) doubles

// Initialize an affine array to arbitrary values
double *set_affine(double a[],double sx,double rx,double ry,double sy,
double tx,double ty)
{
a[0] = sx;
a[1] = rx;
a[2] = ry;
a[3] = sy;
a[4] = tx;
a[5] = ty;
return(a);
}
// Set the affine array to translate by (x,y)
double *set_translate_affine(double t[],double x,double y)
{
t[0] = 1;
t[1] = 0;
t[2] = 0;
t[3] = 1;
t[4] = x;
t[5] = y;
return(t);
}
// Set the affine array to scale the image by sx,sy
double *set_scale_affine(double s[],double sx,double sy)
{
s[0] = sx;
s[1] = 0;
s[2] = 0;
s[3] = sy;
s[4] = 0;
s[5] = 0;
return(s);
}
// Set the affine array to rotate image by 'degrees' clockwise
double *set_rotate_affine(double r[],double degrees)
{
r[0] = cos(DegreesToRadians(fmod(degrees,360.0)));
r[1] = sin(DegreesToRadians(fmod(degrees,360.0)));
r[2] = -sin(DegreesToRadians(fmod(degrees,360.0)));
r[3] = cos(DegreesToRadians(fmod(degrees,360.0)));
r[4] = 0;
r[5] = 0;
return(r);
}
// Multiply two affine arrays and return the result.
double *affine_multiply(double c[],double a[],double b[])
{
c[0] = a[0]*b[0] + a[1]*b[2];
c[1] = a[0]*b[1] + a[1]*b[3];
c[2] = a[2]*b[0] + a[3]*b[2];
c[3] = a[2]*b[1] + a[3]*b[3];
c[4] = a[4]*b[0] + a[5]*b[2] + b[4];
c[5] = a[4]*b[1] + a[5]*b[3] + b[5];
return(c);
}

void test_wand(void) {
// Remember that these operations are done with respect to the
// origin of the image which is the TOP LEFT CORNER.

// Example 1.
// Rotate logo: by 90 degrees (about the origin), scale by 50 percent and
// then move the image 240 in the x direction
{
double s[6],r[6],t[6],temp[6],result[6];
MagickWand *mw = NULL;

MagickWandGenesis();
mw=NewMagickWand();
MagickReadImage(mw,"logo:");

// Set up the affine matrices
// rotate 90 degrees clockwise
set_rotate_affine(r,90);
// scale by .5 in x and y
set_scale_affine(s,.5,.5);
// translate to (240,0)
set_translate_affine(t,240,0);
// now multiply them - note the order in
// which they are specified - in particular beware that
// temp = r*s is NOT necessarily the same as temp = s*r

//first do the rotation and scaling
// temp = r*s
affine_multiply(temp,r,s);
// now the translation
// result = temp*t;
affine_multiply(result,temp,t);

// and then apply the result to the image
MagickDistortImage(mw,AffineProjectionDistortion,
6,(const double *)&result,MagickFalse);

MagickWriteImage(mw,"logo_affine_1.jpg");
if(mw)mw = DestroyMagickWand(mw);
MagickWandTerminus();
}

// Example 2
// Rotate logo: 30 degrees around the point (300,100)
// Since rotation is done around the origin, we must translate
// the point (300,100) up to the origin, do the rotation, and
// then translate back again
//
{
double t1[6],r[6],t2[6],temp[6],result[6];
MagickWand *mw = NULL;

MagickWandGenesis();
mw=NewMagickWand();
MagickReadImage(mw,"logo:");

// Initialize the required affines
// translate (300,100) to origin
set_translate_affine(t1,-300,-100);
// rotate clockwise by 30 degrees
set_rotate_affine(r,30);
// translate back again
set_translate_affine(t2,300,100);

// Now multiply the affine sequence
// temp = t1*r
affine_multiply(temp,t1,r);
// result = temp*t2;
affine_multiply(result,temp,t2);

MagickDistortImage(mw,AffineProjectionDistortion,
6,result,MagickFalse);

MagickWriteImage(mw,"logo_affine_2.jpg");
if(mw)mw = DestroyMagickWand(mw);
MagickWandTerminus();
}

// Example 3
// Reflect the image about a line passing through the origin.
// If the line makes an angle of D degrees with the horizontal
// then this can be done by rotating the image by -D degrees so
// that the line is now (in effect) the x axis, reflect the image
// across the x axis, and then rotate everything back again.
// In this example, rather than just picking an arbitrary angle,
// let's say that we want the "logo:" image to be reflected across
// it's own major diagonal. Although we know the logo: image is
// 640x480 let's also generalize the code slightly so that it
// will still work if the name of the input image is changed.
// If the image has a width "w" and height "h", then the angle between
// the x-axis and the major diagonal is atan(h/w) (NOTE that this
// result is in RADIANS!)
// For this example I will also retain the original dimensions of the
// image so that anything that is reflected outside the borders of the
// input image is lost
{
double r1[6],reflect[6],r2[6],temp[6],result[6];
double w,h,angle_degrees;
MagickWand *mw = NULL;

MagickWandGenesis();
mw=NewMagickWand();
MagickReadImage(mw,"logo:");
w = MagickGetImageWidth(mw);
h = MagickGetImageHeight(mw);

// Just convert the radians to degrees. This way I don't have
// to write a function which sets up an affine rotation for an
// argument specified in radians rather than degrees.
// You can always change this.
angle_degrees = RadiansToDegrees(atan(h/w));
// Initialize the required affines
// Rotate diagonal to the x axis
set_rotate_affine(r1,-angle_degrees);
// Reflection affine (about x-axis)
// In this case there isn't a specific function to set the
// affine array (like there is for rotation and scaling)
// so use the function which sets an arbitrary affine
set_affine(reflect,1,0,0,-1,0,0);
// rotate image back again
set_rotate_affine(r2,angle_degrees);

// temp = r1*reflect
affine_multiply(temp,r1,reflect);
// result = temp*r2;
affine_multiply(result,temp,r2);

MagickDistortImage(mw,AffineProjectionDistortion,
6,result,MagickFalse);

MagickWriteImage(mw,"logo_affine_3.jpg");
if(mw)mw = DestroyMagickWand(mw);
MagickWandTerminus();
}

// Example 4
// Create a rotated gradient
// See: http://www.imagemagick.org/discourse-server/viewtopic.php?f=1&t=12707
// The affine in this one is essentially the same as the one in Example 2 but
// this example has a different use for the result
{
double t1[6],r[6],t2[6],temp[6],result[6];
MagickWand *mw = NULL;
long w,h,nw,nh,gw,gh;
double theta,rad_theta;
char info[256];

// Dimensions of the final rectangle
w = 600;
h = 100;
// angle of clockwise rotation
theta = 15;	// degrees
// Convert theta to radians
rad_theta = DegreesToRadians(theta);
// Compute the dimensions of the rectangular gradient
gw = (long)(w*cos(rad_theta) + h*sin(rad_theta) +0.5);
gh = (long)(w*sin(rad_theta) + h*cos(rad_theta) +0.5);
// Don't let the rotation make the gradient rectangle any smaller
// than the required output
if(gw < w)gw = w;
if(gh < h)gh = h;

MagickWandGenesis();
mw=NewMagickWand();
MagickSetSize(mw,gw,gh);
MagickReadImage(mw,"gradient:white-black");

// Initialize the required affines
// translate centre of gradient to origin
set_translate_affine(t1,-gw/2,-gh/2);
// rotate clockwise by theta degrees
set_rotate_affine(r,theta);
// translate back again
set_translate_affine(t2,gw/2,gh/2);

// Now multiply the affine sequences
// temp = t1*r
affine_multiply(temp,t1,r);
// result = temp*t2;
affine_multiply(result,temp,t2);

MagickDistortImage(mw,AffineProjectionDistortion,6,result,MagickFalse);
// Get the size of the distorted image and crop out the middle
nw = MagickGetImageWidth(mw);
nh = MagickGetImageHeight(mw);
MagickCropImage(mw,w,h,(nw-w)/2,(nh-h)/2);
MagickWriteImage(mw,"rotgrad_2.png");
if(mw)mw = DestroyMagickWand(mw);
MagickWandTerminus();
}

}