Rotation algorithm not producing expected results

Question

I am trying to port some code from IDL to C, and I find myself having to replicate the ROT function. The goal is to rotate a 1024x1024 array of unsigned shorts by an angle in degrees. For the purposes of this project, the angle is very small, less than one degree. The function uses bilinear interpolation.

I tried a backwards approach. For each pixel in the output array, I did a reverse rotation to figure out what coordinate in the input array it would belong to, then used interpolation to figure out what that value would be. I wasn't sure how to go about doing bilinear interpolation if the input grid was skewed; every example of it I've seen assumes that it's orthogonal.

For the rotation, I referred to this:

x' = x * cos(a) + y * sin(a)
y' = y * cos(a) - x * sin(a)

from this: Image scaling and rotating in C/C++

And for the interpolation, I referred to this: http://en.wikipedia.org/wiki/Bilinear_interpolation

Anyway, here's my code:

#define DTOR 0.0174532925

void rotatearray(unsigned short int *inarray, unsigned short int *outarray, int xsize,
                 int ysize, double angle)
{
  //allocate temparray, set to 0                                                        
  unsigned short int *temparray;
  temparray = calloc(xsize*ysize, sizeof(unsigned short int));

  int xf, yf;
  int xi1, xi2, yi1, yi2;
  double xi, yi;
  double x, y;
  double minusangle = (360 - angle)*DTOR;
  unsigned short int v11, v12, v21, v22;

  int goodpixels=0;
  int badpixels=0;

  for(yf=0;yf<ysize;yf++)
    {
      for(xf=0;xf<xsize;xf++)
        {
          //what point in the input grid would map to this output pixel?                
          //(inverse of rotation)                                                       
          xi = (xf+0.5)*cos(minusangle) + (yf+0.5)*sin(minusangle);
          yi = (yf+0.5)*cos(minusangle) - (xf+0.5)*sin(minusangle);

          //Is it within bounds?                                                        
          if ((xi>(0+0.5))&&(xi<xsize-0.5)&&
              (yi>(0+0.5))&&(yi<ysize-0.5))
            {
              //what are the indices of the bounding input pixels?                  
              xi1 = (int)(xi - 0.5);                          
              xi2 = (int)(xi + 0.5);                        
              yi1 = (int)(yi - 0.5);
              yi2 = (int)(yi + 0.5);

              //What position is (x,y) in the bound unit square?
              x = xi - xi1;           
              y = yi - yi1;

              //what are the values of the bounding input pixels?                  
              v11 = inarray[yi1*xsize + xi1];//What are the values of                   
              v12 = inarray[yi2*xsize + xi1];//the bounding input pixels?               
              v21 = inarray[yi1*xsize + xi2];
              v22 = inarray[yi2*xsize + xi2];

              //Do bilinear interpolation                                               
              temparray[yf*xsize + xf] = (unsigned short int)
                (v11*(1-x)*(1-y) + v21*x*(1-y) + v12*(1-x)*y + v22*x*y);
              goodpixels++;
            }
          else{temparray[yf*xsize + xf]=0; badpixels++;}
        }
    }

  //copy to outarray                                                                    
  for(yf=0;yf<ysize;yf++)
    {
      for(xf=0;xf<xsize;xf++)
        {
          outarray[yf*xsize + xf] = temparray[yf*xsize+xf];
        }
    }

  free(temparray);
  return;
}

I tested it by printing several dozen numbers, and comparing it to the same ind of the IDL code, and the results are not at all the same. I'm not sure what more information I can give on that, as I'm not currently able to produce a working image of the array. Do you see any errors in my implementation? Is my reasoning behind the algorithm sound?

EDIT: Here are some selected numbers from the input array; they are identical in the C and IDL programs. What's printed is the x index, followed by the y index, followed by the value at that point.

   0       0      24.0000
 256       0      17.0000
 512       0      23.0000
 768       0      21.0000
1023       0      0.00000
   0     256      19.0000
 256     256      459.000
 512     256      379.000
 768     256      191.000
1023     256      0.00000
   0     512      447.000
 256     512      388.000
 512     512      231.000
 768     512      231.000
1023     512      0.00000
   0     768      286.000
 256     768      378.000
 512     768      249.000
 768     768      205.000
1023     768      0.00000
   0    1023      6.00000
 256    1023      10.0000
 512    1023      11.0000
 768    1023      12.0000
1023    1023      0.00000

This is what the IDL program outputs after rotation:

   0       0      31.0000
 256       0      20.4179
 512       0      20.3183
 768       0      20.0000
1023       0      0.00000
   0     256      63.0000
 256     256      457.689
 512     256      392.406
 768     256      354.140
1023     256      0.00000
   0     512      511.116
 256     512      402.241
 512     512      230.939
 768     512      240.861
1023     512      0.00000
   0     768      296.826
 256     768      377.217
 512     768      218.039
 768     768      277.194
1023     768      0.00000
   0    1023      14.0000
 256    1023      8.00000
 512    1023      9.34906
 768    1023      23.7820
1023    1023      0.00000

And here is the data after rotation using my function:

[0,0]: 0
[256,0]: 44
[512,0]: 276
[768,0]: 299
[1023,0]: 0
[0,256]: 0
[256,256]: 461
[512,256]: 439
[768,256]: 253
[1023,256]: 0
[0,512]: 0
[256,512]: 377
[512,512]: 262
[768,512]: 379
[1023,512]: 0
[0,768]: 0
[256,768]: 340
[512,768]: 340
[768,768]: 198
[1023,768]: 18
[0,1023]: 0
[256,1023]: 0
[512,1023]: 0
[768,1023]: 0
[1023,1023]: 0

I didn't see an immediately useful pattern emerging here to indicate what's going on, which is why I didn't originally include it.

EDIT EDIT EDIT: I believe my mind just suddenly stumbled across the problem! I noticed that the 0,0 pixel never seems to change, and that the 1023,1023 pixel changes the most. Of course this means the algorithm is designed to rotate around the origin, while I'm assuming that the function I seek to imitate is designed to rotate around the center of the image. The answer is still not the same but it is much closer. All I did was change the line

      xi = (xf+0.5)*cos(minusangle) + (yf+0.5)*sin(minusangle);
      yi = (yf+0.5)*cos(minusangle) - (xf+0.5)*sin(minusangle);

to

      xi = (xf-512+0.5)*cos(minusangle) + (yf-512+0.5)*sin(minusangle) + 512;
      yi = (yf-512+0.5)*cos(minusangle) + (xf-512+0.5)*sin(minusangle) + 512;

I'm guessing `DTOR` is degrees to radians, but it might be good to show that value. — macduff, Jul 02 '13 at 20:10
The only information you've provided regarding what's wrong is "...the results are not at all the same". Can you clarify what the error or difference is that you're seeing comparing expected versus seen results? — lurker, Jul 02 '13 at 20:14
I thought it might be information overload, but I guess that's silly. I am adding a list of numbers I printed after doing the rotation in IDL and doing the rotation in C. — Frank Harris, Jul 02 '13 at 20:15
Your C output seems to be `int`, not `double`. Can you double check your format string? — jxh, Jul 02 '13 at 20:23
@jxh, the input data is an array of unsigned shorts, and so I set the output data as such as well. I cast the value to an unsigned short when I did the interpolation. As such, I'm willing to accept a discrepancy of less than 1, but at the moment the discrepancies are much larger than 1. — Frank Harris, Jul 02 '13 at 20:32
They're generated from .RAW images. It's complicated and involves some bitshifting. Furthermore, I printed these same indices in the arrays in my C code and in my IDL code throughout several steps. I found that they were in perfect agreement right up until the rotation. Do you still want the code? — Frank Harris, Jul 02 '13 at 20:47
No, that's not what I meant. Basically, I want to know exactly how to create an array of input that will generate the same input table that you listed, and confirm that I get the same output that you listed. Can you provide the C code that generated the output? — jxh, Jul 02 '13 at 21:38
Debug ideas: A) after `y = yi - yi1;` assert x & y are in 0 to 1 range. B) Fill image with known values (0, 1, or 1/2 max ushort) and look for rotation pattern. C) Verify rotation of 0 works. BTW: should not `x = xi - xi1;` be something like `x = (xi - xi1) / 1024`? Good luck — chux - Reinstate Monica, Jul 02 '13 at 22:13
@chux A) was a VERY useful tip! It totally revealed what a, although not the only, problem was. With the last bit, it's not (xi - xi1)/1024, but xi - xi1 - 0.5. x is supposed to be the difference between the target point and the point of the center of the pixel to the left. That 0.5 corrects for that. I found that rotation of 0 works except for the fact that I check to make that it's not in the outer half-pixel border. I can modify that. I'll try the known value idea next. Thanks for your ideas! — Frank Harris, Jul 03 '13 at 14:46
@jxh I'm not sure how to do that except by sharing some of my group's data, and I'm not sure that we're allowed to do that. I don't think there's a problem with me posting my own code, but the data is another issue. Sorry. — Frank Harris, Jul 03 '13 at 14:51
Does [this](http://stackoverflow.com/questions/620745/c-rotating-a-vector-around-a-certain-point) help perhaps? — meaning-matters, Jul 03 '13 at 17:40
@meaning-matters We couldn't have had better timing. I came to that realization just a moment before I saw this comment. I adjusted my question to reflect the fix. It's still not perfect, but one HUGE error has been fixed! — Frank Harris, Jul 03 '13 at 17:52

Rotation algorithm not producing expected results

0 Answers0