Changeset 41164 for trunk/Ohana/src/opihi/cmd.data/rotate.c

Timestamp:

Nov 27, 2019, 11:22:02 AM (7 years ago)

Author:

eugene

Message:

extend cut function to apply alternate stats to the output vector; extend extract to allow source region to fall beyond bounds of input image; add commands antialias, kapamemory, match1d, mgaussdev, vsigmoid (fit a sigmoid), virls (irls on vector), vwtmean (weighted mean), vtransitions; extend medimage to report variance; fix error in rotate; update some tests; add vgauss -apply option

File:

• trunk/Ohana/src/opihi/cmd.data/rotate.c (modified) (4 diffs)

trunk/Ohana/src/opihi/cmd.data/rotate.c

@@ -3 +3 @@
 int rotate (int argc, char **argv) {
   
-  int i, j, NX, NY, X, Y, Lx, Ly, N;
+  int i, j, NX, NY;
   float *in_buff, *out_buff, *c;
-  double angle, CosAngle, SinAngle, Xo, Yo, dX, dY, fx, fy, x, y, X1, Y1;
+  double Xo, Yo, x, y, X1, Y1;
   double pc11, pc12, pc21, pc22, PC11, PC12, PC21, PC22;
   Buffer *buf;
 
-  Xo = 0;
-  Yo = 0;
-  if ((N = get_argument (argc, argv, "-center"))) {
-    remove_argument (N, &argc, argv);
-    Xo  = atof(argv[N]);
-    remove_argument (N, &argc, argv);
-    Yo  = atof(argv[N]);
-    remove_argument (N, &argc, argv);
-  }
+//  Xo = 0;
+//  Yo = 0;
+//  if ((N = get_argument (argc, argv, "-center"))) {
+//    remove_argument (N, &argc, argv);
+//    Xo  = atof(argv[N]);
+//    remove_argument (N, &argc, argv);
+//    Yo  = atof(argv[N]);
+//    remove_argument (N, &argc, argv);
+//  }
 
   if (argc != 3) {
@@ -221 +221 @@
   }
 
-  angle = atof (argv[2]);
-  CosAngle = cos (angle*RAD_DEG);
-  SinAngle = sin (angle*RAD_DEG);
+  double angle = atof (argv[2]);
+  double CosAngle = cos (angle*RAD_DEG);
+  double SinAngle = sin (angle*RAD_DEG);
   
   gprint (GP_ERR, "rotating: %f %f %f\n", angle, CosAngle, SinAngle);
 
-  Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle);
-  Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle);
-  dX = MAX(0,NY*SinAngle);
-  dY = MAX(0,-NX*SinAngle);
-  /*
-  gprint (GP_ERR, "%f %f  -->  ", Xo, Yo);
-  X1 = Xo*CosAngle - Yo*SinAngle + dX;
-  Y1 =  Xo*SinAngle + Yo*CosAngle + dY;
-  gprint (GP_ERR, "%f %f\n", X1, Y1);
-  */
+  // we are rotating about the center pixel, (NX/2, NY/2),
+  // but we are then putting the result in a new image
+  // of size (Lx,Ly).
+
+  // (x,y) = (i - Nx/2),(j - Ny/2)
+  // (x',y') = R(theta) (x,y)
+  // (I,J) = (x' + Lx/2),(y' + Ly/2)
+
+  int Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle);
+  int Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle);
 
   /* fix reference pixel */
   gfits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo);
   gfits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo);
-  /*
-  X1 = (Xo - dX)*CosAngle - (Yo - dY)*SinAngle;
-  Y1 = (dX - Xo)*SinAngle + (Yo - dY)*CosAngle;
-  */
-  X1 = Xo*CosAngle - Yo*SinAngle + dX;
-  Y1 =  Xo*SinAngle + Yo*CosAngle + dY;
+  X1 = (Xo - NX/2)*CosAngle + (Yo - NY/2)*SinAngle + Lx/2;
+  Y1 = (NX/2 - Xo)*SinAngle + (Yo - NY/2)*CosAngle + Ly/2;
   gfits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
   gfits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);
@@ -255 +251 @@
   gfits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
   gfits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
-  PC11 = pc11*CosAngle - pc21*SinAngle;
-  PC12 = pc12*CosAngle - pc22*SinAngle;
-  PC21 = pc21*CosAngle + pc11*SinAngle;
-  PC22 = pc22*CosAngle + pc12*SinAngle;
+  PC11 = pc11*CosAngle + pc21*SinAngle;
+  PC12 = pc12*CosAngle + pc22*SinAngle;
+  PC21 = pc21*CosAngle - pc11*SinAngle;
+  PC22 = pc22*CosAngle - pc12*SinAngle;
   gfits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
   gfits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
@@ -270 +266 @@
   gfits_create_matrix (&buf[0].header, &buf[0].matrix);
   gfits_print_alt (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
+
   out_buff = (float *)buf[0].matrix.buffer;
   for (j = 0; j < Ly; j++) {
     for (i = 0; i < Lx; i++, out_buff++) {
 
-      x = (i - dX)*CosAngle + (j - dY)*SinAngle;
-      y = (dX - i)*SinAngle + (j - dY)*CosAngle;
-      X = (int) x;
-      Y = (int) y;
-
-      if (X < 0) continue;
-      if (X >= NX - 1) continue;
-      if (Y < 0) continue;
-      if (Y >= NY - 1) continue;
-
-      c = &in_buff[X + NX*Y];
-      fx = x - X;
-      fy = y - Y;
+      float xo = (i - Lx/2);
+      float yo = (j - Ly/2);
+
+      x =  xo*CosAngle + yo*SinAngle;
+      y = -xo*SinAngle + yo*CosAngle;
+
+      int I = x + NX/2;
+      int J = y + NY/2;
+
+      if (I < 0) continue;
+      if (I >= NX - 1) continue;
+      if (J < 0) continue;
+      if (J >= NY - 1) continue;
+
+      c = &in_buff[I + NX*J];
+
+      int X = (int) x;
+      int Y = (int) y;
+      double fx = x - X;
+      double fy = y - Y;
+
       *out_buff = (c[0]*(1-fx) + c[1]*fx)*(1-fy) + (c[NX+1]*fx + c[NX]*(1-fx))*fy;
+//    *out_buff = c[0];
     }
   }