00001 #include "mex.h" /* must appear first */ 00002 #include <stdio.h> 00003 #include <math.h> 00004 #include "mexhead.h" 00005 #include "Doc/rotate_docstring.h" /* autogenerated from this file */ 00006 #include <strings.h> 00007 #include <string.h> // strtok 00008 #include <sys/time.h> // for timing 00009 00010 00011 /************************************************************** 00012 00013 %rotate: rotates the sun for a given number of hours (new version) 00014 % 00015 % y=rotate(x,geom1,geom2,hours,mode,bb); 00016 % * Given a solar image x with pose given by geom1, rotate it backward 00017 % in time by hours hours, into the pose given by geom2. Off-disk and 00018 % unknown values are indicated by NaN. 00019 % * The differential rotation coefficients are defined by this routine. 00020 % The rotation speed depends on latitude, according to: 00021 % omega = coefs(1) + coefs(2) * S + coefs(3) * S^2 00022 % where S = [sin(latitude)]^2, and omega is the apparent 00023 % angular velocity in degrees/day. The equator rotates fastest, 00024 % but this is not assumed in the code. 00025 % * We use obs_v = 0.9865 deg/day and coefs = [14.6,-2.2,0] deg/day. 00026 % * The required mode parameter selects two things, sesw or transposed 00027 % pixel ordering, and also full vs. sparse rotation. These work 00028 % as follows. 00029 % * The normal HMI (and normal MDI) pixel ordering starts in the 00030 % southeast corner, and the first scan line of pixels runs toward the 00031 % southwest corner. This is `sesw' ordering. The transposed ordering 00032 % is `sene'; this ordering is what we used for the JPL MDI processing. 00033 % In particular, the older rotate routine uses the transposed ordering. 00034 % This is implemented via internal stride parameters. 00035 % * The regular full-image rotation looks at every pixel in the 00036 % target image and interpolates a value for each. Or, a 'sparse' 00037 % algorithm can be specified when the rotated input x is mostly zero, 00038 % as for example indicators of active regions. This method does 00039 % a pre-scan of x to determine where the nonzero values are, and 00040 % only computes y in a box bounding the one where x will move to, 00041 % not over the whole disk. It thus saves computation depending on 00042 % the extent of the nonzeros in x (more compact is better). 00043 % * If the bounding box of nonzeros in x is known, it can be given 00044 % as the last argument and that region will be used to restrict the 00045 % pre-scan of x. 00046 % * When mode contains `sparse', uncomputed values are left as zero 00047 % outside this box. Such zero values can be known, on-disk values 00048 % that happen to be zero, or off-disk value, or values that would be 00049 % unknown. In short, no distinction is made between 0 and NaN as 00050 % opposed to when mode contains `full'. But in particular, all 00051 % non-zero, non-NaN values will match up. 00052 % 00053 % Inputs: 00054 % real x(m,n); 00055 % real geom1(5); -- [center_x center_y r_sun beta p0] 00056 % opt real geom2(5) = geom1; -- (empty defaults to above) 00057 % real hours; 00058 % char mode = 'full,sesw' -- required, in: {full,sparse}x{sesw,sene} 00059 % opt int bb = [1 1 m n] -- matlab coordinates 00060 % 00061 % Outputs: 00062 % real y(m,n); 00063 % 00064 % See also: 00065 00066 % implemented as a mex file 00067 % turmon rewrote (again) Sept 2010 00068 ****************************************************************/ 00069 00070 /* constants and globals used for error checking */ 00071 00072 #define NARGIN_MIN 4 /* min number of inputs */ 00073 #define NARGIN_MAX 6 /* max number of inputs */ 00074 #define NARGOUT_MIN 1 /* min number of outputs */ 00075 #define NARGOUT_MAX 1 /* max number of outputs */ 00076 00077 #define ARG_X 0 00078 #define ARG_GEOM1 1 00079 #define ARG_GEOM2 2 00080 #define ARG_HRS 3 00081 #define ARG_MODE 4 00082 #define ARG_BB 5 00083 00084 #define ARG_Y 0 00085 00086 static const char *progname = "rotate"; 00087 #define PROGNAME rotate 00088 static const char *in_specs[NARGIN_MAX] = { 00089 "RM", 00090 "RV(5)", 00091 "RS|RV(5)", // RS permits [], which we want 00092 "RS", 00093 "SV", 00094 "RS|RV(4)"}; // permits [] 00095 static const char *in_names[NARGIN_MAX] = { 00096 "x", 00097 "geom1", 00098 "geom2", 00099 "hours", 00100 "mode", 00101 "bb"}; 00102 static const char *out_names[NARGOUT_MAX] = { 00103 "y"}; 00104 00105 static int verbose = 0; 00106 00107 /* 00108 * Solar rotation constants 00109 * all in deg/day 00110 * first has sidereal rate, less Earth's rotational rate, to make synodic rate 00111 * if you're sensitive to these values, you're probably doing something wrong... 00112 */ 00113 00114 #define ROTATE_CONST_K0 (14.6-0.9865) 00115 #define ROTATE_CONST_K1 (-2.2) 00116 #define ROTATE_CONST_K2 0.0 00117 00118 /* 00119 * extracts various mode parameters from a free-form string 00120 * 00121 * returns 0 if the string was OK, an explanatory tag otherwise 00122 */ 00123 static 00124 char * 00125 extract_mode(char *mode, int *Mfull, int *Msesw) 00126 { 00127 char *word; 00128 char *sep = ", \t"; 00129 00130 if (!mode) return "null"; /* should not happen */ 00131 for (word = strtok(mode, sep); word; word = strtok(NULL, sep)) { 00132 /* full/sparse */ 00133 if (strcasecmp(word, "full" ) == 0) *Mfull = 1; 00134 else if (strcasecmp(word, "sparse") == 0) *Mfull = 0; 00135 /* sesw, or transposed */ 00136 else if (strcasecmp(word, "sesw" ) == 0) *Msesw = 1; 00137 else if (strcasecmp(word, "sene" ) == 0) *Msesw = 0; 00138 /* unrecognized word: not OK */ 00139 else return word; 00140 } 00141 return 0; /* OK */ 00142 } 00143 00144 00145 00146 /*** 00147 *** 00148 *** BASIC ROTATION SUBROUTINES 00149 *** 00150 ***/ 00151 00152 00153 /* 00154 * interpolate in one dimension, returning 00155 * l * i0 + (1-l) * i1 00156 * if possible. If both are NaN, NaN is returned. 00157 * If only i0 is NaN, i1 is returned, unless l = 1. 00158 * The latter indicates no contribution from i1 00159 * is allowed, so NaN must be returned. 00160 * Similarly if only i1 is NaN. 00161 */ 00162 inline 00163 double 00164 interp1(double l, 00165 double i0, 00166 double i1, 00167 double nan) 00168 00169 { 00170 switch ((isnan(i1) ? 2 : 0) + (isnan(i0) ? 1 : 0)) { 00171 case 0: 00172 return l*i0 + (1.0-l)*i1; /* interpolate */ 00173 break; 00174 case 1: 00175 return (l < 1.0) ? i1 : nan; /* i0 is nan, i1 isn't */ 00176 break; 00177 case 2: 00178 return (l > 0.0) ? i0 : nan; /* i1 is nan, i0 isn't */ 00179 break; 00180 case 3: 00181 return nan; /* hopeless */ 00182 break; 00183 } 00184 } 00185 00186 00187 /* 00188 * is this a valid index into an image: 00189 * 0 <= i1 < M1 and 0 <= i2 < M2 00190 * used in do_rotate and find_extent 00191 * 00192 * undef'd at the end of this file 00193 */ 00194 #define ISVALID(i1,i2,M1,M2) (((i1)>=0) && ((i1)<M1) && \ 00195 ((i2)>=0) && ((i2)<M2)) 00196 00197 00198 /* 00199 * do_rotate: perform the rotation, b = rotate(a) 00200 * 00201 * Algorithm is to scan "target" or b image, and at each site (x,y), try 00202 * to map it back to the "source" or a image at (xp,yp). This can't be 00203 * done if (x,y) is not on-disk in the target (in which case, put NaN). 00204 * It also fails when (xp, yp) is not visible on the source image, which 00205 * happens when it has rotated around to the back side. Again, put NaN. 00206 * 00207 * The computation is done only in a bounding box around some region 00208 * of interest in the target image, given by bb[01][xy]. Ordinarily 00209 * the box encompasses the whole image. The sparse mode works by 00210 * narrowing the range of the box. 00211 * 00212 * The mapping (x,y) -> (xp,yp) can be done many ways. A general 00213 * and foolproof way is to associate a site with a position on 00214 * the sphere in three dimensions, say 00215 * p1 = (p1x,p1y,p1z). 00216 * Here, (x,y) is the image plane, and z projects back to the 00217 * observer (z>=0 for the visible side). The map from target 00218 * "backwards" to source means three rotations: 00219 * first, by beta in the y-z plane, then by an amount delta where 00220 * delta = f_kappa(hours,sin(lat)^2), 00221 * kappa = [k0 k1 k2], 00222 * in the x-z plane, and then back by beta in the y-z plane. 00223 * 00224 * In vector notation, this can be written 00225 * p2 = T_{-beta} S_{delta} T_{beta} p1 00226 * where T and S are 3x3 rotation matrices parameterized by the 00227 * given angles. Note that the reverse transformation is 00228 * p1 = T_{-beta} S_{-delta} T_{beta} p2 00229 * since inv(T_{beta}) = T_{-beta}. 00230 * 00231 * In particular, since T_beta leaves x unchanged, 00232 * 1 0 0 00233 * T_beta = 0 cosB sinB 00234 * 0 -sinB cosB 00235 * and similarly for S_delta. Note that delta is a function of 00236 * latitude(p1) == latitude(p2), and that the sine of this latitude 00237 * is just the y component of T_{beta} p1. 00238 * 00239 * Anyway, the algorithm to find (xp,yp) just implements these 00240 * three matrix multiplies, with some calculation in between to 00241 * determine delta. If z<0, the desired point in the target is 00242 * off-disk, so fill in NaN. Otherwise, interpolate the source 00243 * image using bilinear interpolation around the "quadrant" 00244 * [floor(xp),ceil(xp)] X [floor(yp),ceil(yp)]. 00245 */ 00246 00247 static 00248 void 00249 do_rotate(int bb0x, // "bottom" bounding box corner 00250 int bb0y, 00251 int bb1x, // "top" bounding box corner 00252 int bb1y, 00253 int maxx, // image size 00254 int maxy, 00255 int strx, // stride of a and b along x, units of doubles 00256 int stry, // stride along y 00257 double *a, // input image 00258 double *b, // output image 00259 // geom1 00260 double x01, // x center, in C coords, origin at 0 00261 double y01, // y center, in C coords, origin at 0 00262 double R1, // radius, in pixels 00263 double B1, // tilt angle, radians 00264 double P1, // twist angle, radians 00265 // geom2 00266 double x02, // x center, in C coords, origin at 0 00267 double y02, // y center, in C coords, origin at 0 00268 double R2, // radius, in pixels 00269 double B2, // tilt angle, radians 00270 double P2, // twist angle, radians 00271 // velocity coeffs 00272 double k0, // order-0 coefficient, radians 00273 double k1, // order-1 coefficient, radians 00274 double k2) // order-2 coefficient, radians 00275 00276 { 00277 int x, y, off0; // target ("b") image coordinates, offset 00278 double p1x, p1y, p1z; // point P1 coords in 3D 00279 double p2x, p2y, p2z; // point P2 coords in 3D 00280 double temp; // workspace slot 00281 double cosD, sinD; // rotation parameters 00282 int xs0, ys0; // position of "lower" pixel in source (a) image 00283 double a00, a01, a10, a11; // aXY = src image @ x=xs0+X, yt=ys0+Y 00284 double lambdaX, lambdaY; // wts: in 1D, b(xt) = lambda*a0 + (1-lambda)*a1 00285 const double nan = mxt_getnand(); // cache nan 00286 // cache geom1 00287 const double cosB1 = cos(B1); 00288 const double sinB1 = sin(B1); 00289 const double cosP1 = cos(P1); 00290 const double sinP1 = sin(P1); 00291 const double R1_R2 = R1/R2; 00292 // cache geom2 00293 const double cosB2 = cos(B2); 00294 const double sinB2 = sin(B2); 00295 const double cosP2 = cos(P2); 00296 const double sinP2 = sin(P2); 00297 const double R2_inv = 1.0/R2; 00298 struct timeval tv1, tv2; 00299 00300 gettimeofday(&tv1, NULL); 00301 /* loop over all target image pixels, filling each one in */ 00302 for (x = bb0x; x <= bb1x; x++) 00303 for (y = bb0y; y <= bb1y; y++) { 00304 /* find target image ("2") vector coordinates P2 = (px2, py2, pz2) */ 00305 p2x = x - x02; 00306 p2y = y - y02; 00307 /* computing: temp = R*R - p2x*p2x - p2y*p2y, but want to avoid 00308 * subtracting large squares, e.g. when p2x ~= +/- R, which is not 00309 * numerically good. We can do so for one subtraction (only) via 00310 * a*a - b*b = (a - b)*(a + b). 00311 * That way, one factor or the other will be close to zero, before 00312 * squaring. 00313 * We can only rewrite using x or y, but not both. So we should 00314 * check which is closer to +/- R, and expand that term. This costs 00315 * too much time, so for now just rewrite the x term only, 00316 * which improves matters half the time. 00317 */ 00318 temp = (R2 + p2x) * (R2 - p2x) - p2y*p2y; /* R*R-x*x = (R-x)*(R+x) */ 00319 if (temp < 0) { b[x*strx+y*stry] = nan; continue; } /* off-disk */ 00320 p2z = sqrt(temp); /* the positive sqrt is the visible half-sphere */ 00321 /* rotate P2 into P1 in the xy plane (the p-angle) -- z is unchanged */ 00322 p1x = cosP2 * p2x + sinP2 * p2y; 00323 p1y = -sinP2 * p2x + cosP2 * p2y; 00324 p1z = p2z; 00325 /* rotate by B2 in yz plane (x axis fixed): P2 = rot(P1) */ 00326 p2x = p1x; /* no change to x */ 00327 p2y = cosB2 * p1y + sinB2 * p1z; 00328 p2z = -sinB2 * p1y + cosB2 * p1z; 00329 /* find the rotation angle due to the time difference */ 00330 temp = R2_inv * p2y; // extract sin(lat) from p2y 00331 temp *= temp; // temp is now sin(lat)^2 00332 temp = k0 + temp*(k1 + k2*temp); // k0 + k1*temp + k2*temp^2 00333 /* temp is now the rotation angle delta, find its sine and cosine */ 00334 cosD = cos(temp); 00335 sinD = sin(temp); 00336 /* rotate by delta in xz plane (y fixed): P1 = rot(P2) */ 00337 p1x = cosD * p2x + sinD * p2z; 00338 p1y = p2y; /* no change to y */ 00339 p1z = -sinD * p2x + cosD * p2z; 00340 /* de-rotate by B1 in yz plane (x fixed): P2 = rot(P1) */ 00341 p2x = p1x; /* no change to x */ 00342 p2y = cosB1 * p1y - sinB1 * p1z; 00343 p2z = sinB1 * p1y + cosB1 * p1z; 00344 /* check for points absent from the front side of the source (a) image */ 00345 if (p2z < 0) { b[x*strx+y*stry] = nan; continue; } /* off-disk */ 00346 /* de-rotate in xy plane by P1, scale by R1/R2, translate back to (x01, y01) */ 00347 p1x = R1_R2 * (cosP1 * p2x - sinP1 * p2y) + x01; 00348 p1y = R1_R2 * (sinP1 * p2x + cosP1 * p2y) + y01; 00349 /* in essence b(x,y) = a(p1x,p1y). But must approximate 00350 a(p1x,p1y), so use bilinear interpolation on the 00351 four pixels surrounding this point, ignoring NaN's */ 00352 /* "lower" corner of box in source image containing remapped (x,y) */ 00353 xs0 = floor(p1x); 00354 ys0 = floor(p1y); 00355 /* Find the box of enclosing pixels; if off-image, substitute a NaN */ 00356 off0 = xs0 * strx + ys0 * stry; // offset from a(0,0) to a(xs0,ys0) 00357 a00 = ISVALID(xs0, ys0, maxx, maxy) ? a[off0 ] : nan; 00358 a01 = ISVALID(xs0, ys0+1, maxx, maxy) ? a[off0+ stry] : nan; 00359 a10 = ISVALID(xs0+1, ys0, maxx, maxy) ? a[off0+strx ] : nan; 00360 a11 = ISVALID(xs0+1, ys0+1, maxx, maxy) ? a[off0+strx+stry] : nan; 00361 /* find the weights associated with smaller-indexed pixels */ 00362 lambdaX = 1 - (p1x - xs0); 00363 lambdaY = 1 - (p1y - ys0); 00364 /* bilinear interpolation within the box to find the new pixel */ 00365 b[x*strx+y*stry] = interp1(lambdaX, 00366 interp1(lambdaY, a00, a01, nan), 00367 interp1(lambdaY, a10, a11, nan), nan); 00368 } /* end for (x,y) */ 00369 gettimeofday(&tv2, NULL); 00370 int dtval = (((long) tv2.tv_sec) - ((long) tv1.tv_sec)) * 1000000 + 00371 (((int) tv2.tv_usec) - ((int) tv1.tv_usec)); 00372 if (verbose) 00373 mexPrintf(" finishing do_rotate: dt = %.1f ms\n", dtval/1000.0); 00374 } 00375 00376 /*** 00377 *** 00378 *** SPARSE BOX-FINDING SUBROUTINES 00379 *** 00380 *** The remaining routines work together to find one (m,n) bounding 00381 *** box, within the target, that contains the rotated version of 00382 *** all non-zero, non-nan, on-disk source pixels. 00383 ***/ 00384 00385 /* 00386 * update_lon: update longitude components of a lat-lon bounding 00387 * box by passing the previous interval through a latitude-dependent 00388 * mapping. 00389 * Background: It is convenient to use interval arithmetic and note that 00390 * continuous functions of intervals in R yield intervals in R (e.g., Rudin, 00391 * theorem 4.22). 00392 * So, the interval we start with will map to another interval, and we 00393 * don't have to do anything fancy along the way, except to remember 00394 * the current interval. (This throws away some constraints that, if 00395 * exploited, could further constrain the interval -- but this will always 00396 * supply an upper bound on interval location.) 00397 * The mapping is 00398 * [lonp0,lonp1] = [lon0,lon1] + [del0,del1] 00399 * where [del0,del1] = f_kappa(sin^2([lat0,lat1])) 00400 * Above f_kappa is the quadratic, parameterized by three coeffs k0, k1, k2: 00401 * f_kappa(a) = k0 + a k1 + a*a k2 00402 * where (above) a = sin^2(lat). 00403 * Since [lat0,lat1] is an interval, so is sin^2 of it, and so is f_kappa 00404 * of that interval, and indeed so is the overall mapping to get lonp. 00405 * The only tricky point is that all functions aren't monotonic, so we need 00406 * to be careful about (for example) 00407 * sin^2([-pi/4,+pi/4]) != [sin^2(-pi/4),sin^2(pi/4)] 00408 * In general this could be hard, but not for quadratics. In particular, if 00409 * g(.) is quadratic, then the endpoints of the interval g([a,b]) are 00410 * A = g(a) or B = g(b), or C = g(c), where c is the point in [a,b] (if any) 00411 * that extremizes the quadratic. So we just need to take min(A,B,C) and 00412 * max(A,B,C) to find the two endpoints of the interval, where we understand 00413 * that we set C = B (say) if c is not in [a,b]. 00414 */ 00415 /* note, these evaluate their arguments several times */ 00416 /* min of three arguments */ 00417 #define min3(a,b,c) ( ((a) <= (b)) ? ( ((a) <= (c)) ? (a) : (c) ) : \ 00418 ( ((b) <= (c)) ? (b) : (c) ) ) 00419 /* max of three arguments */ 00420 #define max3(a,b,c) ( ((a) >= (b)) ? ( ((a) >= (c)) ? (a) : (c) ) : \ 00421 ( ((b) >= (c)) ? (b) : (c) ) ) 00422 00423 static 00424 void 00425 update_lon(double *lon0p, 00426 double *lon1p, 00427 double k0, 00428 double k1, 00429 double k2, 00430 double lat0, 00431 double lat1, 00432 double lon0, 00433 double lon1) 00434 { 00435 double ptA, ptB, ptC; /* holding places for possible endpoints */ 00436 double alpha0, alpha1; /* for the range of sin^2(.) mapping */ 00437 double delta0, delta1; /* for the range of f(sin^2(.)) mapping */ 00438 double alphaT; /* temporary alpha value */ 00439 00440 /* first mapping: [alpha0,alpha1] = sin^2([lat0,lat1]) */ 00441 /* (note, sin(.) is easy because it's monotonic for latitude angles) */ 00442 ptA = sin(lat0); ptA *= ptA; /* one endpoint of the interval */ 00443 ptB = sin(lat1); ptB *= ptB; /* other endpoint */ 00444 ptC = (lat0 * lat1 < 0) ? 0.0 : ptB; /* ptC = sin^2(0) = 0 if 0 in [lat0,lat1] */ 00445 alpha0 = min3(ptA,ptB,ptC); 00446 alpha1 = max3(ptA,ptB,ptC); 00447 /* second mapping: [delta0,delta1] = f_kappa([alpha0,alpha1]) */ 00448 ptA = -(k0 + k1 * alpha0 + k2 * alpha0 * alpha0); 00449 ptB = -(k0 + k1 * alpha1 + k2 * alpha1 * alpha1); 00450 ptC = ptB; /* if unchanged below, this will be a no-op */ 00451 if (k2 != 0) { 00452 alphaT = -k1/(2 * k2); /* the extremizing (max or min) alpha */ 00453 /* test if alphaT in [alpha0,alpha1] */ 00454 if ((alpha0 - alphaT) * (alpha1 - alphaT) < 0) 00455 ptC = -(k0 + k1 * alphaT + k2 * alphaT * alphaT); /* ptC = f_kappa(alphaT) */ 00456 } 00457 delta0 = min3(ptA,ptB,ptC); 00458 delta1 = max3(ptA,ptB,ptC); 00459 /* finally, [lon0p,lon1p] = [lon0,lon1] + [delta0,delta1] */ 00460 *lon0p = lon0 + delta0; 00461 *lon1p = lon1 + delta1; 00462 } 00463 00464 /* don't leave these macros hanging around */ 00465 #undef min3 00466 #undef max3 00467 00468 00469 /* 00470 * pointmap: map a point in (lat,lon) coordinates on a sphere, at tilt 00471 * angle B given by (sinB, cosB), to (x,y,z) coordinates. 00472 */ 00473 static 00474 void 00475 pointmap(double *px, double *py, double *pz, /* output: point coords */ 00476 double s_lat, double c_lat, /* lat sine and cosine */ 00477 double s_lon, double c_lon, /* lon sine and cosine */ 00478 double sinB, double cosB, 00479 double sinP, double cosP) 00480 { 00481 // intermediate point coords 00482 double qx, qy, qz; 00483 double rx, ry, rz; 00484 00485 qx = c_lat * s_lon; 00486 qy = s_lat; 00487 qz = c_lat * c_lon; 00488 /* rotate by -beta in yz plane (x axis fixed): R = rot(Q) */ 00489 rx = qx; /* no change to x */ 00490 ry = cosB * qy - sinB * qz; 00491 rz = sinB * qy + cosB * qz; 00492 /* rotate by -p in xy plane (z axis fixed): P = rot(R) */ 00493 *px = cosP * rx - sinP * ry; 00494 *py = sinP * rx + cosP * ry; 00495 *pz = rz; // no change to z 00496 } 00497 00498 /* 00499 * pointunmap: map a point in (x,y,z) coordinates on a sphere, at tilt 00500 * angle B given by (sinB, cosB), to (lat,lon) coordinates. 00501 */ 00502 static 00503 void 00504 pointunmap(double *lat, double *lon, /* output: point coords */ 00505 double px, double py, double pz, 00506 double sinB, double cosB, 00507 double sinP, double cosP) 00508 { 00509 // intermediate point coords 00510 double qx, qy, qz; 00511 double rx, ry, rz; 00512 00513 /* rotate by +p in xy plane (z axis fixed): R = rot(P) */ 00514 rx = cosP * px + sinP * py; 00515 ry = -sinP * px + cosP * py; 00516 rz = pz; // no change to z 00517 /* rotate by +beta in yz plane (x axis fixed): Q = rot(R) */ 00518 qx = rx; /* no change to x */ 00519 qy = cosB * ry + sinB * rz; 00520 qz = -sinB * ry + cosB * rz; 00521 /* find lat and lon */ 00522 *lon = atan2(qx, qz); /* p2=(1,0,0), at limb, has lon=pi/2 */ 00523 *lat = asin(qy); 00524 } 00525 00526 /* 00527 * find_enclosing_xy_bb: find a bounding box in x-y coordinates, say 00528 * [bb0x,bb1x] X [bb0y,bb1y] 00529 * that encloses a bounding box in lat-lon coordinates given by 00530 * [lat0,lat1] X [lon0, lon1]. 00531 * The lat-lon coordinates are on a sphere at tilt angle beta given 00532 * by (sinB, cosB). The lon coordinates given as input do not have 00533 * to be in (-pi,pi). This is because they arise from adding an 00534 * angular displacement to another pair of lon coordinates, and either 00535 * or both results may land outside (-pi,pi). Because we want to preserve 00536 * the ordering (lon0 < lon1), it is important to not wrap them. (And 00537 * there is no compelling reason to do so.) 00538 * 00539 * There may be an analytic way to do this. This code uses that we can 00540 * easily traverse the edges of the lat-lon patch, and that the x-y 00541 * bounding box will need to be only as big as area enclosed by the 00542 * boundary. It therefore divides the boundary into four sides, and 00543 * sets up a loop to march along these sides in "nstep" steps. At each 00544 * step, the four points along the four boundaries are computed, and 00545 * their (x,y,z) coords are found from the (lat,lon). The overall span 00546 * of x and y is retained as the bounding box. 00547 * 00548 * There are two wrinkles, both related to edge effects. 00549 * 00550 * First, the traversal is discrete and may be interrupted by the limb. 00551 * If it is, we may not go out far enough onto the limb, in (lat,lon) 00552 * coordinates, to encompass the final on-boundary pixel. Hence, 00553 * we indicate by "expand" whether the disk boundary was crossed. If 00554 * it is, the entire boundary is expanded by one pixel. (Discretization 00555 * can cause other problems, but they are more of a second-order nature 00556 * than this one -- e.g., a curved boundary may be traversed in steps, 00557 * but the fact that it pokes slightly above a pixel border not detected, 00558 * because it happens between steps.) 00559 * 00560 * Second, wrapping issues may make the x/y extrema not occur on the 00561 * boundary. Take beta = 0 for simplicity. If a square patch at disk 00562 * center is rotated to the limb, it may extend partway around the limb 00563 * at the equator. 00564 * 00565 * In this case, both traversed (longitudinal) boundaries will 00566 * curve back around the limb at values of x strictly less than 00567 * unity. In particular, the boundary does not include the point 00568 * where the equator crosses behind the limb, but the patch does. 00569 * This point has x=+/-1 and is thus the true extreme of the patch 00570 * in the x direction. 00571 * 00572 * What has happened? The extrema of x and y can occur either at 00573 * the boundary of the patch (which is the "expected case," and will 00574 * be found in the loop), or at the unconstrained extrema: 00575 * (x,y,z) = (+/-1, 0, 0) (for x) 00576 * (x,y,z) = ( 0, +/-1, 0) (for y) 00577 * These unconstrained extrema may be in the patch, but not on the boundary, 00578 * so we have to check them. (There are no other special cases.) 00579 * 00580 * To take care of these cases, we explicitly compute the (lat,lon) coordinates 00581 * of these four unconstrained extrema, checking if they fall into the lat-lon 00582 * bounding box. If they do, then the lat/lon bounding box is adjusted 00583 * accordingly. 00584 */ 00585 00586 // macro updates the 4 BB sides based on a point (px,py,pz) 00587 #define UPDATE_BB \ 00588 if (pts) {pts[inx++] = px; pts[inx++] = py; pts[inx++] = pz;} \ 00589 if (pz >= 0) { \ 00590 if (bb0x > px) bb0x = px; \ 00591 if (bb1x < px) bb1x = px; \ 00592 if (bb0y > py) bb0y = py; \ 00593 if (bb1y < py) bb1y = py; \ 00594 } else { \ 00595 expand = 1; \ 00596 } 00597 00598 00599 static 00600 void 00601 find_enclosing_xy_bb(int *P_bb0x, // BB (x,y) boundary 00602 int *P_bb0y, 00603 int *P_bb1x, 00604 int *P_bb1y, 00605 double lat0, // BB (lat,lon) boundary... 00606 double lat1, // lat in -pi/2 .. pi/2 00607 double lon0, // lon, in radians, unrestricted 00608 double lon1, 00609 int maxx, // image domain 00610 int maxy, 00611 double x0, // x center, in C coords, origin at 0 00612 double y0, // y center, in C coords, origin at 0 00613 double R, // radius, in pixels 00614 double B, // tilt angle, radians 00615 double P) // twist angle, radians 00616 { 00617 // counters and dimensions 00618 int n; // count position on enclosing boundary 00619 int expand; // expand bb by 1 if it wraps to the far side 00620 double degree_per_pixel; 00621 int npixel; 00622 double d_lat; // step size along boundaries 00623 double d_lon; 00624 // point parameters 00625 double px, py, pz; // their (x,y,z) versions 00626 double lat, lon; // varying positions along boundaries 00627 double s_lat, c_lat, s_lon, c_lon; // sin,cos of lat,lon (varying with n) 00628 double bb0x, bb0y, bb1x, bb1y; // abbreviate the bb boundaries 00629 // cache geom 00630 const double cosB = cos(B); 00631 const double sinB = sin(B); 00632 const double cosP = cos(P); 00633 const double sinP = sin(P); 00634 // cache lat/lon bb 00635 const double s_lat0 = sin(lat0); /* sin, cos of lat0, lon0 */ 00636 const double c_lat0 = cos(lat0); 00637 const double s_lon0 = sin(lon0); 00638 const double c_lon0 = cos(lon0); 00639 const double s_lat1 = sin(lat1); /* sin, cos of lat1, lon1 */ 00640 const double c_lat1 = cos(lat1); 00641 const double s_lon1 = sin(lon1); 00642 const double c_lon1 = cos(lon1); 00643 // allow export of points 00644 double *pts = NULL; // pointer to triples: (x, y, z) 00645 int inx; 00646 int npts; 00647 00648 // ensure there is at least one test point per image pixel along each edge 00649 // there are R pixels per radian at disk center (2*pi*R pixels / 2*pi radians) 00650 // so, number of pixels is (pixels/radian) * radians = R * radians 00651 npixel = ceil(R * max(lat1-lat0, lon1-lon0)); 00652 if (npixel > max(maxx, maxy)) npixel = max(maxx, maxy); 00653 if (npixel < 10) npixel = 10; 00654 d_lat = (lat1 - lat0) / npixel; 00655 d_lon = (lon1 - lon0) / npixel; 00656 npts = (npixel + 1) * 4; // npixel+1 points X 4 sides 00657 pts = mxCalloc(npts*3, sizeof(double)); // each point has 3 coords 00658 // assume we did not cross onto the back side 00659 expand = 0; 00660 // right up until the end of this function, bb is on the unit sphere 00661 bb0x = bb0y = 1.0; bb1x = bb1y = -1.0; // min <- highval; max <- lowval 00662 /* move points 1,2,3,4 simultaneously over all four boundaries */ 00663 for (inx=n=0, lat=lat0, lon=lon0; n<=npixel; n++, lat+=d_lat, lon+=d_lon) { 00664 /* sine, cosine of lat and lon markers (used by all points) */ 00665 s_lat = sin(lat); c_lat = cos(lat); 00666 s_lon = sin(lon); c_lon = cos(lon); 00667 /* from lat-lon to xyz: point1 = (lat0, lon) */ 00668 pointmap(&px, &py, &pz, s_lat0, c_lat0, s_lon, c_lon, sinB, cosB, sinP, cosP); 00669 UPDATE_BB; 00670 /* from lat-lon to xyz: point2 = (lat1, lon) */ 00671 pointmap(&px, &py, &pz, s_lat1, c_lat1, s_lon, c_lon, sinB, cosB, sinP, cosP); 00672 UPDATE_BB; 00673 /* from lat-lon to xyz: point3 = (lat, lon0) */ 00674 pointmap(&px, &py, &pz, s_lat, c_lat, s_lon0, c_lon0, sinB, cosB, sinP, cosP); 00675 UPDATE_BB; 00676 /* from lat-lon to xyz: point4 = (lat, lon1) */ 00677 pointmap(&px, &py, &pz, s_lat, c_lat, s_lon1, c_lon1, sinB, cosB, sinP, cosP); 00678 UPDATE_BB; 00679 } 00680 /* Adjust bounding box to allow for disk-edge (z=0) boundary. 00681 * Note that the computed lon of these key points is in -pi..pi, so it 00682 * must be adjusted to fall in to the [lon0,lon1] interval, which may 00683 * lie outside -pi..pi. */ 00684 /* (1,0,0) */ 00685 pointunmap(&lat, &lon, 1.0, 0.0, 0.0, sinB, cosB, sinP, cosP); 00686 while (lon < lon0) lon += 2*M_PI; /* wrap lon to best fit in [lon0,lon1] */ 00687 while (lon > lon1) lon -= 2*M_PI; 00688 if (lat0 <= lat && lat <= lat1 && lon0 <= lon && lon <= lon1) 00689 bb1x = 1.0; 00690 /* (-1,0,0) */ 00691 pointunmap(&lat, &lon, -1.0, 0.0, 0.0, sinB, cosB, sinP, cosP); 00692 while (lon < lon0) lon += 2*M_PI; /* wrap lon to best fit in [lon0,lon1] */ 00693 while (lon > lon1) lon -= 2*M_PI; 00694 if (lat0 <= lat && lat <= lat1 && lon0 <= lon && lon <= lon1) 00695 bb0x = -1.0; 00696 /* (0,1,0) */ 00697 pointunmap(&lat, &lon, 0.0, 1.0, 0.0, sinB, cosB, sinP, cosP); 00698 while (lon < lon0) lon += 2*M_PI; /* wrap lon to best fit in [lon0,lon1] */ 00699 while (lon > lon1) lon -= 2*M_PI; 00700 if (lat0 <= lat && lat <= lat1 && lon0 <= lon && lon <= lon1) 00701 bb1y = 1.0; 00702 /* (0,-1,0) */ 00703 pointunmap(&lat, &lon, 0.0, -1.0, 0.0, sinB, cosB, sinP, cosP); 00704 while (lon < lon0) lon += 2*M_PI; /* wrap lon to best fit in [lon0,lon1] */ 00705 while (lon > lon1) lon -= 2*M_PI; 00706 if (lat0 <= lat && lat <= lat1 && lon0 <= lon && lon <= lon1) 00707 bb0y = -1.0; 00708 00709 /* define the BB corners in full-disk (x,y) coords... */ 00710 bb0x = floor(bb0x * R + x0) - expand; 00711 bb0y = floor(bb0y * R + y0) - expand; 00712 bb1x = ceil(bb1x * R + x0) + expand; 00713 bb1y = ceil(bb1y * R + y0) + expand; 00714 /* ...and ensure the BB corners are on-image */ 00715 if (bb0x < 0) bb0x = 0; 00716 if (bb0y < 0) bb0y = 0; 00717 if (bb1x >= maxx) bb1x = maxx-1; 00718 if (bb1y >= maxy) bb1y = maxy-1; 00719 // plug in output BB (output BB are int's) 00720 *P_bb0x = (int) bb0x; 00721 *P_bb0y = (int) bb0y; 00722 *P_bb1x = (int) bb1x; 00723 *P_bb1y = (int) bb1y; 00724 00725 if (pts) { 00726 mxt_put_matrix("bdry", -1, pts, 3, npts); 00727 mxFree(pts); 00728 } 00729 00730 } 00731 00732 /* 00733 * Given image a, find the lat-lon bounding box of all its non-zero, non-nan pixels. 00734 * 00735 * We consider only the pixels in a in the rectangle [bb0m,bb1m]x[bb0n,bb1n]. 00736 * We use a workspace "b" which is image-sized, and assumed zero-filled. 00737 */ 00738 static 00739 void 00740 find_lat_lon_bb(double *p_lat0, // low side of (lat,lon) boundary... 00741 double *p_lon0, // lon, in radians, unrestricted 00742 double *p_lat1, // lat in -pi/2 .. pi/2 00743 double *p_lon1, 00744 int bb0x, // bottom corner of BB to search 00745 int bb0y, 00746 int bb1x, // top corner of BB 00747 int bb1y, 00748 int maxx, // input size (# points) 00749 int maxy, 00750 int strx, // stride of a and b along x, units of doubles 00751 int stry, // stride along y 00752 double *a, // input is maxx X maxy 00753 char *b, // workspace is same size, already zero-filled 00754 // geom 00755 double x0, // x center, in C coords, origin at 0 00756 double y0, // y center, in C coords, origin at 0 00757 double R, // radius, in pixels 00758 double B, // tilt angle, radians 00759 double P) // twist angle, radians 00760 { 00761 int x, y; // image coordinates 00762 int dx, dy; // coordinates for neighbors 00763 int off1, off2; // array offsets 00764 double p1x, p1y, p1z; // point P1 coords in 3D 00765 double p2x, p2y, p2z; // point P2 coords in 3D 00766 double temp; // workspace slot 00767 // cache geom 00768 const double cosB = cos(B); 00769 const double sinB = sin(B); 00770 const double cosP = cos(P); 00771 const double sinP = sin(P); 00772 double lat, lon; // current latitude and longitude 00773 double lat0, lon0, lat1, lon1; // lat-lon BB while accumulating 00774 struct timeval tv1, tv2; // crude timings 00775 00776 #if 0 00777 double *zs = mxCalloc(maxx*maxy, sizeof(double)); 00778 double *lats = mxCalloc(maxx*maxy, sizeof(double)); 00779 double *lons = mxCalloc(maxx*maxy, sizeof(double)); 00780 #endif 00781 00782 gettimeofday(&tv1, NULL); 00783 /* Loop over source image pixels, finding a (lat,lon) BB */ 00784 lon0 = 10.0; lon1 = -10.0; /* will be reset on first iter */ 00785 lat0 = R+1; lat1 = -(R+1); /* will be reset on first iter */ 00786 for (x = bb0x; x <= bb1x; x++) 00787 for (y = bb0y; y <= bb1y; y++) { 00788 off1 = x*strx + y*stry; 00789 /* skip zero or NaN (unknown or off-disk) pixels */ 00790 if (a[off1] == 0 || isnan(a[off1])) continue; 00791 /* ensure (x,y) is itself on-disk */ 00792 p1x = x - x0; 00793 p1y = y - y0; 00794 temp = (R + p1x) * (R - p1x) - p1y*p1y; // R*R-x*x = (R-x)*(R+x) 00795 // zs[x*strx + y*stry] = temp; 00796 if (temp < 0) continue; // (x,y) is off-disk: skip 00797 /* found a nonzero, on-disk pixel: get its nbrs lat-lon */ 00798 for (dx = x-1; dx <= x+1; dx++) 00799 for (dy = y-1; dy <= y+1; dy++) { 00800 off2 = dx*strx + dy*stry; 00801 if (!ISVALID(dx,dy,maxx,maxy) || b[off2]) 00802 continue; 00803 else 00804 b[off2] = 1; 00805 /* find original vector coordinates P2 = (px2, py2, pz2) */ 00806 p2x = dx - x0; 00807 p2y = dy - y0; 00808 temp = (R + p2x) * (R - p2x) - p2y*p2y; // R*R-x*x = (R-x)*(R+x) 00809 /* temp is the residual which goes into the p1z component */ 00810 if (temp < 0) { 00811 p2z = -sqrt(-temp); // barely off-disk - send to back side 00812 temp = R / sqrt(p2x*p2x + p2y*p2y + p2z*p2z); // need ||P2|| = 1 00813 p2x *= temp; p2y *= temp; p2z *= temp; // rescale 00814 } else { 00815 p2z = sqrt(temp); /* totally OK -- on visible side */ 00816 } 00817 /* rotate by p-angle in xy plane (z axis fixed): P1 = rot(P2) */ 00818 p1x = cosP * p2x + sinP * p2y; 00819 p1y = -sinP * p2x + cosP * p2y; 00820 p1z = p2z; 00821 /* rotate by +beta in yz plane (x axis fixed): P2 = rot(P1) */ 00822 p2x = p1x; /* no change to x */ 00823 p2y = cosB * p1y + sinB * p1z; 00824 p2z = -sinB * p1y + cosB * p1z; 00825 /* find (lat,lon) of the point */ 00826 lon = atan2(p2x, p2z); // p2=(1,0,0), at limb, has lon=pi/2 00827 lat = p2y; // take the arcsin out of the loop: need range only 00828 00829 // lats[off2] = asin(lat/R); 00830 // lons[off2] = lon; 00831 00832 /* update bounding box */ 00833 if (lat < lat0) lat0 = lat; 00834 if (lat > lat1) lat1 = lat; 00835 if (lon < lon0) lon0 = lon; 00836 if (lon > lon1) lon1 = lon; 00837 } /* for dx, dy */ 00838 } /* for x, y */ 00839 /* map latitude range through (monotone) arcsin function */ 00840 lat0 = asin(lat0 / R); 00841 lat1 = asin(lat1 / R); 00842 /* set up output pointers */ 00843 *p_lat0 = lat0; 00844 *p_lon0 = lon0; 00845 *p_lat1 = lat1; 00846 *p_lon1 = lon1; 00847 gettimeofday(&tv2, NULL); 00848 int dtval = (((long) tv2.tv_sec) - ((long) tv1.tv_sec)) * 1000000 + 00849 (((int) tv2.tv_usec) - ((int) tv1.tv_usec)); 00850 if (verbose) 00851 mexPrintf(" finishing find_lat_lon_bb, dt = %.1f ms\n", dtval/1000.0); 00852 00853 #if 0 00854 if (lats) { 00855 mxt_put_matrix("zs", -1, zs, maxx, maxy); 00856 mxt_put_matrix("lats", -1, lats, maxx, maxy); 00857 mxt_put_matrix("lons", -1, lons, maxx, maxy); 00858 mxFree(zs); 00859 mxFree(lats); 00860 mxFree(lons); 00861 } 00862 #endif 00863 } 00864 00865 00866 /* 00867 * find_extent: finds an (m,n) bounding box around the location, 00868 * in the target image, of all nontrivial source image pixels. 00869 * This works in three steps: 00870 * 1: find a (lat,lon) bounding box around all nontrivial source 00871 * pixels. (Nontrivial == on-disk, non-zero, non-NaN.) The step 00872 * is done in this routine. 00873 * 2: propagate the (lat,lon) bounding box to the target image 00874 * using the rules for differential rotation. The target BB is 00875 * held in (latp,lonp). This step is done primarily by update_lon, 00876 * since latp = lat. 00877 * 3: find a (m,n) bounding box around the (latp,lonp) bounding 00878 * box. This step is done in find_enclosing_xy_bb. 00879 * 00880 * The differential rotation and the discretization make this 00881 * process tricky. Steps 2 and 3 are mostly explained in their 00882 * subroutines. 00883 * As for step 1, the basic principle is to loop over all 00884 * source pixels, skipping trivial pixels. The remainder 00885 * are mapped to (lat,lon) coordinates. The max and min values 00886 * of both are maintained. The lon that is used here is zero at disk 00887 * center, so that it wraps around from pi->-pi not at either limb, 00888 * but at the far back side of the sun. This is done so that 00889 * splotches that are near the limb will not have pixels at lon = pi 00890 * and lon = -pi. Such values would force the lon range (which must 00891 * be held as a simple interval for efficiency and clarity) to 00892 * be -pi..pi. This would not be incorrect, just inefficient. 00893 * 00894 * The main refinement has to do with how interpolation is done by 00895 * do_rotate around NaN values, particularly at the limb. If a 00896 * target image pixel lands at (say) (500.1,61.1) in the source, 00897 * the closest source pixel is at (500,61). Suppose this pixel 00898 * is NaN. The interpolated value will be chosen using (501,62), 00899 * or (500,62), or (501,61), if any or all are non-NaN: 00900 * 00901 * 61 62 00902 * 500 NaN 4 00903 * 501 NaN 3 00904 * 00905 * Thus, from do_rotate's point of view, (500,62) is NOT the "last 00906 * available" nontrivial pixel. Actually, (500,61+epsilon) is, 00907 * where epsilon is arbitrarily close to zero. Looking now at 00908 * the src-> target mapping, this means that every valid source 00909 * pixel (x,y) implies we know a valid value at (x+dx,y+dy) for 00910 * all (dx,dy) in [-1,1] X [-1,1]. 00911 * 00912 * This can be observed in practice. 00913 * 00914 * One solution: Every time a nontrivial pixel is found, loop around 00915 * the 8 neighbors of the pixel, and compute their lat-lon's as well, 00916 * to add them to the bounding box. This would imply a lot more work, 00917 * so we save the "included pixels" into a scratch array b, which is 00918 * assumed ZERO at the outset. So, we don't waste time re-including 00919 * the same pixels. 00920 * 00921 * The wrinkle is that some of these surrogate pixels will fall off-disk, 00922 * when (x+dx,y+dy) is farther than R from (n0,m0) -- even though (x,y) 00923 * was not. This typically means that the square patch corresponding to 00924 * the surrogate pixel falls only partly on the visible side of the disk: 00925 * the center happens to be off-disk. 00926 * 00927 * They can't simply be discarded. An ad hoc solution is to wrap the patch 00928 * around the disk. This is done by setting a negative "z" coordinate 00929 * (which will have a small value relative to R). In this case, it is 00930 * vital that its coordinates in 3D have the right norm, i.e. the 00931 * "wrapped" center must fall on the sphere. (This is obvious.) 00932 * 00933 * The simplest test cases show the value of this approach. Take 00934 * beta = 0 and an active region partly intersecting one limb. Rotate 00935 * it slightly towards disk center. Without this approach, the lon 00936 * part of the bounding box will not go up to pi/2 because no pixel 00937 * center happens to fall right on the disk edge. With it, the edge 00938 * pixels that are only partly visible are wrapped back so their 00939 * centers are on the back side, giving a lon > pi/2. The bounding 00940 * box will then be wide enough to contain all valid pixels. 00941 * 00942 */ 00943 static 00944 void 00945 find_extent(int *bb0x, // bottom corner of BB (in/out) 00946 int *bb0y, 00947 int *bb1x, // top corner of BB (in/out) 00948 int *bb1y, 00949 int maxx, // input size (# points) 00950 int maxy, 00951 int strx, // stride of a and b along x, units of doubles 00952 int stry, // stride along y 00953 double *a, // input 00954 char *b, // workspace arrives zero-filled 00955 // geom1 00956 double x01, // x center, C coords, origin at 0 00957 double y01, // y center, C coords, origin at 0 00958 double R1, // radius, in pixels 00959 double B1, // tilt angle, radians 00960 double P1, // twist angle, radians 00961 // geom2 00962 double x02, // x center, C coords, origin at 0 00963 double y02, // y center, C coords, origin at 0 00964 double R2, // radius, in pixels 00965 double B2, // tilt angle, radians 00966 double P2, // twist angle, radians 00967 // velocity coeffs 00968 double k0, // order-0 coefficient, radians 00969 double k1, // order-1 coefficient, radians 00970 double k2) // order-2 coefficient, radians 00971 00972 { 00973 double lat0, lon0, lat1, lon1; // lat-lon BB at start time 00974 double lat0p, lon0p, lat1p, lon1p; // lat-lon BB at end time 00975 00976 #ifdef DEBUG 00977 printf("FE1: mnBB = (%d,%d) (%d,%d)\n", *bb0x, *bb0y, *bb1x, *bb1y); 00978 #endif 00979 /* Step 1: Find a (lat,lon) BB within the source image `a' */ 00980 find_lat_lon_bb(&lat0, &lon0, &lat1, &lon1, 00981 *bb0x, *bb0y, *bb1x, *bb1y, 00982 maxx, maxy, strx, stry, 00983 a, b, 00984 x01, y01, R1, B1, P1); 00985 #ifdef DEBUG 00986 printf("FE2: llBB = (%f,%f) (%f,%f)\n", 00987 lat0*180/M_PI, lon0*180/M_PI, lat1*180/M_PI, lon1*180/M_PI); 00988 #endif 00989 /* no active pixels were found: empty BB: eliminate the boundary case */ 00990 if (lat0 > lat1 || lon0 > lon1) { 00991 *bb0x = *bb0y = *bb1x = *bb1y = 0; 00992 return; 00993 } 00994 00995 /* Step 2: Convert source lat-lon BB to a target lat-lon BB containing it */ 00996 /* step 2a: latitude restrictions do not change */ 00997 lat0p = lat0; 00998 lat1p = lat1; 00999 /* step 2b: find longitude restrictions (more tricky) */ 01000 update_lon(&lon0p, &lon1p, k0, k1, k2, lat0, lat1, lon0, lon1); 01001 #ifdef DEBUG 01002 printf("FE3: llBB = (%f,%f) (%f,%f)\n", 01003 lat0p*180/M_PI, lon0p*180/M_PI, lat1p*180/M_PI, lon1p*180/M_PI); 01004 printf(" lon0: (%f -> %f); lon1: (%f -> %f)\n del_lon: %f, %f %s\n", 01005 lon0*180/M_PI, lon0p*180/M_PI, 01006 lon1*180/M_PI, lon1p*180/M_PI, 01007 (lon1 - lon0)*180/M_PI, (lon1p - lon0p)*180/M_PI, 01008 ((lon1 - lon0) < (lon1p - lon0p)) ? " " : "!!!"); 01009 #endif 01010 01011 /* Step 3: Identify an enclosing BB in the target in (x,y) coords */ 01012 find_enclosing_xy_bb(bb0x, bb0y, bb1x, bb1y, 01013 lat0p, lat1p, lon0p, lon1p, 01014 maxx, maxy, 01015 x02, y02, R2, B2, P2); 01016 #ifdef DEBUG 01017 printf("FE4: mnBBp = (%d,%d) (%d,%d)\n", *bb0x, *bb0y, *bb1x, *bb1y); 01018 #endif 01019 } 01020 01021 /* gateway routine */ 01022 01023 #ifdef StaticP /* undefined under mex */ 01024 StaticP 01025 #endif 01026 void 01027 mexFunction(int nlhs, 01028 mxArray *plhs[], 01029 int nrhs, 01030 const mxArray *prhs[]) 01031 { 01032 int m, n; // original matlab sizes 01033 int maxx, maxy; // sizes, transposed if needed 01034 int strx, stry; // strides in images along x and y dims 01035 int bb0m, bb0n, bb1m, bb1n; // bounding box coords (original matlab) 01036 int bb0x, bb0y, bb1x, bb1y; // bounding box coords (transposed if needed) 01037 char *mode, *word; // label for converted mode string and substring 01038 int mode_full; // normal = 1, sparse = 0 01039 int mode_sesw; // sesw = 1, transposed = 0 01040 double datamin, datamax; 01041 char errstr[120]; 01042 double dt; // time interval 01043 double k0, k1, k2; // rotation coefficients (see below) 01044 double *src, *targ; // pointers into "x" for 2d indexing 01045 char *scratch; // scratchpad of 0/1 indicators 01046 double *g; // abbreviate head of geometry tuple 01047 double x01, y01, r1, b1, p1; // geom1 01048 double x02, y02, r2, b2, p2; // geom2 01049 struct timeval tv1, tv2; // crude timings 01050 01051 gettimeofday(&tv1, NULL); 01052 /* Hook for introspection (function signature, docstring) */ 01053 if (nrhs < 0) { 01054 plhs[0] = mxt_PackSignature((mxt_Signature) (-nrhs), 01055 NARGIN_MIN, NARGIN_MAX, 01056 NARGOUT_MIN, NARGOUT_MAX, 01057 in_names, in_specs, out_names, docstring); 01058 return; 01059 } 01060 /* check number of arguments */ 01061 if ((nrhs < NARGIN_MIN) || (nrhs > NARGIN_MAX)) 01062 mexErrMsgTxt((snprintf(errstr, sizeof(errstr), 01063 "%s: Expect %d <= input args <= %d", 01064 progname, NARGIN_MIN, NARGIN_MAX), errstr)); 01065 if ((nlhs < NARGOUT_MIN) || (nlhs > NARGOUT_MAX)) 01066 mexErrMsgTxt((snprintf(errstr, sizeof(errstr), 01067 "%s: Expect %d <= output args <= %d", 01068 progname, NARGOUT_MIN, NARGOUT_MAX), errstr)); 01069 mexargparse(nrhs, prhs, in_names, in_specs, NULL, progname); 01070 01071 /* get size of input image x, needed below */ 01072 m = (int) mxGetM(prhs[ARG_X]); 01073 n = (int) mxGetN(prhs[ARG_X]); 01074 01075 /****************************************************************** 01076 * Read and check arguments in order supplied 01077 ******************************************************************/ 01078 01079 /* 01080 * image geometry 01081 */ 01082 // source image (x) geometry 01083 g = mxGetPr(prhs[ARG_GEOM1]); // abbreviate geom1 01084 x01 = g[0] - 1; // x0, in C coordinates 01085 y01 = g[1] - 1; // y0, in C coordinates 01086 r1 = g[2]; // disk radius 01087 b1 = g[3]*M_PI/180.0; // deg -> rad 01088 p1 = g[4]*M_PI/180.0; // deg -> rad 01089 // target image (y) geometry 01090 if (mxGetNumberOfElements(prhs[ARG_GEOM2]) > 0) 01091 g = mxGetPr(prhs[ARG_GEOM2]); // use geom2 instead of geom1 01092 x02 = g[0] - 1; // x0, in C coordinates 01093 y02 = g[1] - 1; // y0, in C coordinates 01094 r2 = g[2]; // disk radius 01095 b2 = g[3]*M_PI/180.0; // deg -> rad 01096 p2 = g[4]*M_PI/180.0; // deg -> rad 01097 01098 /* 01099 * Rotation angle coefficients, a function of t 01100 * 01101 * omitting units, the basic formula is: 01102 * given: solar rotation coefficients k, constant, sin(lat)^2, sin(lat)^4 01103 * observer's velocity (sidereal), obs_v 01104 * time interval, dt 01105 * compute: rotation angle delta 01106 * delta = dt * (k(0) + k(1)*sin(lat)^2 + k(2)*sin(lat)^4 - obs_v) 01107 * = dt * ([k(0)-obs_v] + k(1)*sin(lat)^2 + k(2)*sin(lat)^4) 01108 * We put everything in common units (radians/day) and multiply by dt. 01109 * The resulting coeffs need only to be multiplied by the powers of 01110 * sin(lat) to get the rotation angle delta in radians. 01111 */ 01112 dt = mxGetScalar(prhs[ARG_HRS]) / 24.0; /* time, now in days */ 01113 /* note, k0 below has the sidereal velocity component also, 01114 * all three start in units of degrees/day and are converted, deg/day to rad/day */ 01115 k0 = ROTATE_CONST_K0 * M_PI / 180.0; 01116 k1 = ROTATE_CONST_K1 * M_PI / 180.0; 01117 k2 = ROTATE_CONST_K2 * M_PI / 180.0; 01118 /* multiply through by the time in days to obtain radians */ 01119 k0 *= dt; k1 *= dt; k2 *= dt; 01120 01121 /* 01122 * `mode' input 01123 * (required) 01124 */ 01125 if ((mode = mxArrayToString(prhs[ARG_MODE])) == NULL) 01126 mexErrMsgTxt((snprintf(errstr, sizeof(errstr), 01127 "%s: bad mode (non-string?). " 01128 "Could not convert mode arg to string.", 01129 progname), errstr)); 01130 /* no default values -- user must specify both explicitly */ 01131 mode_full = mode_sesw = -1; 01132 if ((word = extract_mode(mode, &mode_full, &mode_sesw)) != NULL) 01133 mexErrMsgTxt((snprintf(errstr, sizeof(errstr), 01134 "%s: bad word <%s> in mode <%.80s>", 01135 progname, word, /* need a fresh copy of the mode */ 01136 mxArrayToString(prhs[ARG_MODE])), errstr)); 01137 if (mode_full < 0 || mode_sesw < 0) 01138 mexErrMsgTxt((snprintf(errstr, sizeof(errstr), 01139 "%s: need full/sparse, sesw/sene in mode <%.80s>", 01140 progname, /* need a fresh copy of the mode */ 01141 mxArrayToString(prhs[ARG_MODE])), errstr)); 01142 mxFree(mode); // done with the mode string 01143 01144 /* 01145 * `bb' input 01146 * (optional) 01147 */ 01148 if (nrhs <= ARG_BB || mxGetNumberOfElements(prhs[ARG_BB]) == 0) { 01149 // not supplied or 0 elements is OK, just take full box 01150 bb0m = 0; bb0n = 0; bb1m = m-1; bb1n = n-1; // C indexing 01151 } else if (mxGetNumberOfElements(prhs[ARG_BB]) == 4) { 01152 // source bounding box was given explicitly 01153 double *bb = mxGetPr(prhs[ARG_BB]); // abbreviation 01154 bb0m = bb[0]; 01155 bb0n = bb[1]; 01156 bb1m = bb[2]; 01157 bb1n = bb[3]; 01158 if ((bb0m < 1 || bb0m > m) || 01159 (bb1m < 1 || bb1m > m) || 01160 (bb0n < 1 || bb0n > n) || 01161 (bb1n < 1 || bb1n > n) || 01162 (bb0m > bb1m || bb0n > bb1n)) 01163 mexErrMsgTxt((snprintf(errstr, sizeof(errstr), 01164 "%s: Bounding box `bb' was empty or out-of-range", 01165 progname), errstr)); 01166 // convert matlab to C indexing 01167 bb0m -= 1; bb0n -= 1; bb1m -= 1; bb1n -= 1; 01168 } else { 01169 mexErrMsgTxt((snprintf(errstr, sizeof(errstr), 01170 "%s: Expect `bb', if nonempty, to be of length 4", 01171 progname), errstr)); 01172 } 01173 01174 /* create output image, which needs to be zero-filled. 01175 * we do allocation after the argument checking 01176 */ 01177 if ((plhs[ARG_Y] = mxCreateDoubleMatrix(m, n, mxREAL)) == NULL) 01178 mexErrMsgTxt((snprintf(errstr, sizeof(errstr), 01179 "%s: Failed to allocate space for output array", 01180 progname), errstr)); 01181 01182 /* 01183 * Depending on pixel-order mode, set up strides and bounding boxes in 01184 * (x,y) coordinates from their (m,n) settings. 01185 * Note that the (x0,y0) positions, in the geom1 and geom2 inputs, 01186 * correspond to x and to y, so they do not need to be swapped. 01187 */ 01188 if (mode_sesw) { 01189 // x (E-W axis) varies fastest -- x corresponds to matlab m coordinate 01190 strx = 1; stry = m; // x varies fastest, x has m points... 01191 maxx = m; maxy = n; // y has n points 01192 bb0x = bb0m; bb0y = bb0n; 01193 bb1x = bb1m; bb1y = bb1n; 01194 } else { 01195 // y (N-S axis) varies fastest -- x corresponds to matlab n coordinate 01196 strx = m; stry = 1; // y varies fastest, y has m points... 01197 maxx = n; maxy = m; // and x has n points 01198 bb0x = bb0n; bb0y = bb0m; 01199 bb1x = bb1n; bb1y = bb1m; 01200 } 01201 01202 // convenience for argument passing 01203 src = mxGetPr(prhs[ARG_X]); 01204 targ = mxGetPr(plhs[ARG_Y]); 01205 01206 /* Depending on sparse/full mode setting, establish bounding box coords 01207 * in source image (C conventions) */ 01208 if (!mode_full) { 01209 // sparse version 01210 // note, scratchpad must be filled with zeros 01211 if ((scratch = calloc(m*n, sizeof(char))) == NULL) 01212 mexErrMsgTxt((snprintf(errstr, sizeof(errstr), 01213 "%s: Failed to allocate space for scratchpad matrix", 01214 progname), errstr)); 01215 // note, bb* has been set up above as a box within `x' 01216 find_extent(&bb0x, &bb0y, &bb1x, &bb1y, // given as input, received as output 01217 maxx, maxy, // size 01218 strx, stry, // strides in units of doubles 01219 src, scratch, // input, and scratchpad 01220 x01, y01, r1, b1, p1, // geom1 01221 x02, y02, r2, b2, p2, // geom2 01222 k0, k1, k2); // coefficients in radians 01223 free(scratch); 01224 } 01225 01226 /* Rotate the input x to fill in the output y */ 01227 do_rotate(bb0x, bb0y, bb1x, bb1y, // target image bounding box 01228 maxx, maxy, // sizes 01229 strx, stry, // strides in units of doubles 01230 src, targ, 01231 x01, y01, r1, b1, p1, // geom1 01232 x02, y02, r2, b2, p2, // geom2 01233 k0, k1, k2); // coefficients in radians 01234 01235 /* initialize range of output */ 01236 getrange(prhs[ARG_X], &datamin, &datamax); 01237 setrange(plhs[ARG_Y], datamin, datamax); 01238 gettimeofday(&tv2, NULL); 01239 int dtval = (((long) tv2.tv_sec) - ((long) tv1.tv_sec)) * 1000000 + 01240 (((int) tv2.tv_usec) - ((int) tv1.tv_usec)); 01241 if (verbose) 01242 mexPrintf("Exiting %s, dt = %.1f ms\n", progname, dtval/1000.0); 01243 } 01244 01245 01246 /* this macro was defined above */ 01247 #undef ISVALID 01248 01249 01250 /* Hook for generic tail matter */ 01251 #ifdef MEX2C_TAIL_HOOK 01252 #include "mex2c_tail.h" 01253 #endif 01254