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

/*
	Originally inspired by:
https://imagemagick.dev.org.tw/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:
	https://imagemagick.dev.org.tw/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: https://imagemagick.dev.org.tw/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();
}

}