Changeset 5624
- Timestamp:
- Nov 29, 2005, 4:00:11 PM (20 years ago)
- Location:
- trunk/psLib
- Files:
-
- 9 edited
-
src/astro/psCoord.c (modified) (10 diffs)
-
src/astro/psCoord.h (modified) (2 diffs)
-
src/imageops/psImageConvolve.c (modified) (3 diffs)
-
src/imageops/psImageConvolve.h (modified) (2 diffs)
-
src/math/psMinimize.c (modified) (15 diffs)
-
src/math/psPolynomial.c (modified) (8 diffs)
-
src/math/psStats.c (modified) (2 diffs)
-
src/mathtypes/psImage.c (modified) (1 diff)
-
test/sys/verified/tst_psMemory.stderr (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/astro/psCoord.c
r5588 r5624 10 10 * @author GLG, MHPCC 11 11 * 12 * @version $Revision: 1.9 3$ $Name: not supported by cvs2svn $13 * @date $Date: 2005-11- 23 23:54:43$12 * @version $Revision: 1.94 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2005-11-30 02:00:00 $ 14 14 * 15 15 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 61 61 should rename this. 62 62 63 To derive this transformation, start with the simple 2-by-2 matrix inversion 64 based on discriminants. This will invert the 65 (x_2, y_2) = matrix(a b c d) * vector (x, y) 66 where there are no constant terms. Then you substitute 67 x_2 = x_1 - e 68 y_2 = y_1 - f 69 for (x_2, y_2) to get the desired inverse. 63 XXX: Use the ADD version which is based on determinants. 70 64 *****************************************************************************/ 65 66 // XXX EAM : below is the code using the standard matrix representation. 67 // note that this inversion requires x->nX == 1, y->nY == 1 and 68 // x->nY <= 1, y->nX <= 1 71 69 psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform) 72 70 { … … 130 128 // printf("HMMM: out->y: (%f %f %f)\n", out->y->coeff[0][0], out->y->coeff[1][0], out->y->coeff[0][1]); 131 129 130 131 // unless the cross terms are available, set these matrix elements to 0 132 psF64 r12 = 0.0; 133 if (transform->x->nY == 1) { 134 r12 = transform->x->coeff[0][1]; 135 } 136 psF64 r21 = 0.0; 137 if (transform->y->nX == 1) { 138 r21 = transform->y->coeff[1][0]; 139 } 140 psF64 r11 = transform->x->coeff[1][0]; 141 psF64 r22 = transform->y->coeff[0][1]; 142 psF64 xo = transform->x->coeff[0][0]; 143 psF64 yo = transform->y->coeff[0][0]; 144 145 psF64 invDet = 1.0 / (r11 * r22 - r12 * r21); 146 147 // apply the results back to the polynomials 148 out->x->coeff[0][0] = -invDet * (r22 * xo - r12 * yo); 149 out->y->coeff[0][0] = -invDet * (r11 * yo - r21 * xo); 150 out->x->coeff[1][0] = +invDet * r22; 151 out->y->coeff[0][1] = +invDet * r11; 152 if (transform->x->nY == 1) { 153 out->x->coeff[0][1] = -invDet * r12; 154 } 155 if (transform->y->nX == 1) { 156 out->y->coeff[1][0] = -invDet * r21; 157 } 132 158 return(out); 133 159 } … … 327 353 XXX: Private Function. 328 354 329 piNormalize(): take an input angle in radians and convert it to the range 330 0:2*PI. 355 piNormalize(): take an input angle in radians and convert it to the range 0:2*PI. 331 356 *****************************************************************************/ 332 357 psF32 piNormalize(psF32 angle) … … 380 405 PS_ASSERT_PTR_NON_NULL(projection, NULL); 381 406 382 psF64 theta = 0.0; 383 psF64 phi = 0.0; 407 psF64 phi, theta; 408 psF64 sinDp, cosDp, sinAlpha, cosAlpha, sinDelta, cosDelta; 409 psF64 sinTheta, cosPhiCT, sinPhiCT, zeta; 410 411 bool zenithal = (projection->type == PS_PROJ_TAN) ||(projection->type == PS_PROJ_SIN); 384 412 385 413 // Allocate return value … … 391 419 } 392 420 393 // Convert to projection spherical coordinate system 394 theta = asin( sin(coord->d)*sin(projection->D) + 395 cos(coord->d)*cos(projection->D)*cos(coord->r-projection->R)); 396 phi = atan2( -1.0*cos(coord->d)*sin(coord->r-projection->R), 397 sin(coord->d)*cos(projection->D) - cos(coord->d)*sin(projection->D)*cos(coord->r-projection->R) ); 421 if (zenithal) { 422 sinDp = sin(projection->D); 423 cosDp = cos(projection->D); 424 sinAlpha = sin(coord->r-projection->R); 425 cosAlpha = cos(coord->r-projection->R); 426 sinDelta = sin(coord->d); 427 cosDelta = cos(coord->d); 428 429 sinTheta = sinDelta*sinDp + cosDelta*cosDp*cosAlpha; 430 cosPhiCT = sinDelta*cosDp - cosDelta*sinDp*cosAlpha; 431 sinPhiCT = -cosDelta*sinAlpha; 432 } else { 433 phi = coord->r - projection->R; 434 theta = coord->d - projection->D; 435 } 398 436 399 437 // Perform the specified projection 400 // Gnomonic projection 401 if (projection->type == PS_PROJ_TAN) { 402 out->x = (cos(theta)*sin(phi))/sin(theta); 403 out->y = (-1.0*cos(theta)*cos(phi))/sin(theta); 438 switch (projection->type) { 439 case PS_PROJ_TAN: 440 // Gnomonic projection 441 out->x = +sinPhiCT / sinTheta; 442 out->y = -cosPhiCT / sinTheta; 443 break; 444 case PS_PROJ_SIN: 404 445 // Othrographic projection 405 } else if (projection->type == PS_PROJ_SIN) { 406 out->x = cos(theta)*sin(phi); 407 out->y = -1.0*cos(theta)*cos(phi); 446 out->x = +sinPhiCT; 447 out->y = -cosPhiCT; 448 break; 449 case PS_PROJ_AIT: 408 450 // Hammer-Aitoff projection 409 } else if ( projection->type == PS_PROJ_AIT) { 410 psF64 zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0))); 451 zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0))); 411 452 out->x = 2.0*zeta*cos(theta)*sin(phi/2.0); 412 453 out->y = zeta*sin(theta); 454 break; 455 case PS_PROJ_PAR: 413 456 // Parabolic projection 414 } else if ( projection->type == PS_PROJ_PAR) {415 457 out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0); 416 458 out->y = M_PI*sin(theta/3.0); 417 } else {459 default: 418 460 psError(PS_ERR_BAD_PARAMETER_TYPE, true, 419 461 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN, … … 424 466 425 467 // Apply plate scales 426 out->x *= projection->Xs;427 out->y *= projection->Ys;468 out->x /= projection->Xs; 469 out->y /= projection->Ys; 428 470 429 471 // Return output … … 444 486 PS_ASSERT_PTR_NON_NULL(coord, NULL); 445 487 PS_ASSERT_PTR_NON_NULL(projection, NULL); 488 489 psF64 rho = 0.0; 490 psF64 sinTheta = 0.0; 491 psF64 cosTheta = 0.0; 492 psF64 sinPhi = 0.0; 493 psF64 cosPhi = 0.0; 446 494 447 495 psF64 theta = 0.0; … … 458 506 459 507 // Remove plate scales 460 // XXX: Verify this. EAM suggested we do a multiply, however that does 461 // not make sense if we also do the multiply in psProject(). 462 psF64 x = coord->x/projection->Xs; 463 psF64 y = coord->y/projection->Ys; 508 psF64 x = coord->x*projection->Xs; 509 psF64 y = coord->y*projection->Ys; 510 psF64 R = sqrt(x*x + y*y); 511 512 bool zenithal = (projection->type == PS_PROJ_TAN) ||(projection->type == PS_PROJ_SIN); 464 513 465 514 // Perform inverse projection 466 // Gnonomic deprojection 467 if ( projection->type == PS_PROJ_TAN) { 468 phi = atan(-1.0*x/y); 469 theta = atan(1.0/sqrt(x*x+y*y)); 515 switch (projection->type) { 516 case PS_PROJ_TAN: 517 // Gnonomic deprojection 518 rho = sqrt (1 + R*R); 519 sinTheta = 1 / rho; 520 cosTheta = R / rho; 521 sinPhi = (R == 0) ? 0.0 : +x / R; 522 cosPhi = (R == 0) ? 1.0 : -y / R; 523 break; 524 case PS_PROJ_SIN: 470 525 // Orhtographic deprojection 471 } else if ( projection->type == PS_PROJ_SIN) { 472 phi = atan((-1.0*x)/y); 473 theta = atan( sqrt(1.0-(x*x+y*y)) / sqrt(x*x+y*y)); 526 sinTheta = sqrt (1 - R*R); 527 cosTheta = R; 528 sinPhi = (R == 0) ? 0.0 : +x / R; 529 cosPhi = (R == 0) ? 1.0 : -y / R; 530 break; 531 case PS_PROJ_AIT: 474 532 // Hammer-Aitoff deprojection 475 } else if ( projection->type == PS_PROJ_AIT) { 476 psF64 z = sqrt(1.0 - ((x/4.0)*(x/4.0)) - ((y/2.0)*(y/2.0))); 477 phi = 2.0*atan((z*x) / (2.0*(2.0*z*z-1.0)) ); 478 theta = asin(y*z); 533 // XXX EAM : need range check on z^2 : must be > 0 534 // XXX EAM : old code, ADD, and elixir code are discrepant re x/4, y/2 535 rho = sqrt(1.0 - PS_SQR(x/4.0) - PS_SQR(y/2.0)); 536 phi = 2.0*atan2((2.0*rho*rho-1.0), x*rho); 537 theta = asin(y*rho); 538 break; 539 case PS_PROJ_PAR: 479 540 // Parabolic deprojection 480 } else if ( projection->type == PS_PROJ_PAR) { 481 psF64 rho = y/M_PI; 541 rho = y/M_PI; 482 542 phi = x/(1.0 - 4.0*rho*rho); 483 543 theta = 3.0*asin(rho); 484 // Invalid deprojection type485 } else {544 break; 545 default: 486 546 psError(PS_ERR_BAD_PARAMETER_TYPE, true, 487 547 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN, … … 491 551 } 492 552 493 // Convert from projection spherical coordinates 494 out->d = asin( sin(theta)*sin(projection->D) + 495 cos(theta)*cos(projection->D)*cos(phi) ); 496 out->r = projection->R + atan2( -1.0*cos(theta)*sin(phi), 497 sin(theta)*cos(projection->D) - 498 cos(theta)*sin(projection->D)*cos(phi) ); 553 if (zenithal) { 554 psF64 sinDp = sin(projection->D); 555 psF64 cosDp = cos(projection->D); 556 557 // Convert from projection spherical coordinates 558 psF64 delta = asin(sinTheta*sinDp + cosTheta*cosDp*cosPhi); 559 psF64 sinAlphaF = -cosTheta*sinPhi; 560 psF64 cosAlphaF = -cosTheta*cosPhi*sinDp + sinTheta*cosDp; 561 562 out->d = delta; 563 out->r = atan2(sinAlphaF, cosAlphaF) + projection->R; 564 } else { 565 out->r = phi + projection->R; 566 out->d = theta + projection->D; 567 } 499 568 500 569 // Return sphere coordinate -
trunk/psLib/src/astro/psCoord.h
r5542 r5624 1 1 /** @file psCoord.h 2 *3 * @brief Contains basic coordinate transformation definitions and operations4 *5 * This file defines the basic types for astronomical coordinate6 * transformation7 *8 * @ingroup CoordinateTransform9 *10 *11 * @author GLG, MHPCC12 *13 * @version $Revision: 1.47$ $Name: not supported by cvs2svn $14 * @date $Date: 2005-11-18 19:39:29$15 *16 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii17 */2 * 3 * @brief Contains basic coordinate transformation definitions and operations 4 * 5 * This file defines the basic types for astronomical coordinate 6 * transformation 7 * 8 * @ingroup CoordinateTransform 9 * 10 * 11 * @author GLG, MHPCC 12 * 13 * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ 14 * @date $Date: 2005-11-30 02:00:00 $ 15 * 16 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii 17 */ 18 18 19 19 #ifndef PS_COORD_H … … 323 323 psPlaneTransform *transform ///< transform to invert 324 324 ); 325 psPlaneTransform *p_psPlaneTransformLinearInvert_MHPCC( 326 psPlaneTransform *transform ///< transform to invert 327 ); 325 328 326 329 -
trunk/psLib/src/imageops/psImageConvolve.c
r5573 r5624 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1.2 7$ $Name: not supported by cvs2svn $8 * @date $Date: 2005-11- 22 20:15:35$7 * @version $Revision: 1.28 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2005-11-30 02:00:07 $ 9 9 * 10 10 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 487 487 } 488 488 489 voidpsImageSmooth (psImage *image,489 bool psImageSmooth (psImage *image, 490 490 double sigma, 491 491 double Nsigma) … … 569 569 PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, 570 570 typeStr); 571 } 572 } 571 return false; 572 } 573 } 574 return true; 573 575 } 574 576 -
trunk/psLib/src/imageops/psImageConvolve.h
r4920 r5624 7 7 * @author Robert DeSonia, MHPCC 8 8 * 9 * @version $Revision: 1.1 2$ $Name: not supported by cvs2svn $10 * @date $Date: 2005- 08-31 02:07:11$9 * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-11-30 02:00:07 $ 11 11 * 12 12 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 141 141 * Applies a circularly symmetric Gaussian smoothing first in x and then in y 142 142 * directions with just a vector. This process is 2N faster than 2D convolutions (in general). 143 * 144 * @return bool TRUE if successful, otherwise FALSE 143 145 */ 144 voidpsImageSmooth(146 bool psImageSmooth( 145 147 psImage *image, ///< the image to be smoothed 146 148 double sigma, ///< the width of the smoothing kernel in pixels -
trunk/psLib/src/math/psMinimize.c
r5530 r5624 10 10 * @author EAM, IfA 11 11 * 12 * @version $Revision: 1.14 5$ $Name: not supported by cvs2svn $13 * @date $Date: 2005-11- 16 23:06:19 $12 * @version $Revision: 1.146 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2005-11-30 02:00:09 $ 14 14 * 15 15 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 1441 1441 } 1442 1442 1443 // here is the definition for BuildSums4D. ifdef'ed away until it is used 1444 // by psPolynomial4DFit.. 1445 # if (0) 1446 /****************************************************************************** 1447 BuildSums4D(sums, x, y, z, t, nXterm, nYterm, nZterm, nTterm). equiv to 1448 BuildSums2D(). The result is returned as a double **** 1449 *****************************************************************************/ 1450 static double ****BuildSums4D( 1451 psF64 ****sums, 1452 psF64 x, 1453 psF64 y, 1454 psF64 z, 1455 psF64 t, 1456 psS32 nXterm, 1457 psS32 nYterm, 1458 psS32 nZterm, 1459 psS32 nTterm) 1460 { 1461 psS32 nXsum = 0; 1462 psS32 nYsum = 0; 1463 psS32 nZsum = 0; 1464 psS32 nTsum = 0; 1465 psF64 xSum = 1.0; 1466 psF64 ySum = 1.0; 1467 psF64 zSum = 1.0; 1468 psF64 tSum = 1.0; 1469 1470 nXsum = 2*nXterm; 1471 nYsum = 2*nYterm; 1472 nZsum = 2*nZterm; 1473 nTsum = 2*nTterm; 1474 if (sums == NULL) { 1475 sums = (psF64 ****) psAlloc (nXsum*sizeof(psF64)); 1476 for (int i = 0; i < nXsum; i++) { 1477 sums[i] = (psF64 ***) psAlloc (nYsum*sizeof(psF64)); 1478 for (int j = 0; j < nYsum; j++) { 1479 sums[i][j] = (psF64 **) psAlloc (nZsum*sizeof(psF64)); 1480 for (int k = 0; k < nZsum; k++) { 1481 sums[i][j][k] = (psF64 *) psAlloc (nTsum*sizeof(psF64)); 1482 } 1483 } 1484 } 1485 } 1486 // careful with this function: there is no size checking and realloc for reuse 1487 1488 tSum = 1.0; 1489 for (int m = 0; m < nTsum; m++) { 1490 zSum = tSum; 1491 for (int k = 0; k < nZsum; k++) { 1492 ySum = zSum; 1493 for (int j = 0; j < nYsum; j++) { 1494 xSum = ySum; 1495 for (int i = 0; i < nXsum; i++) { 1496 sums[i][j][k][m] = xSum; 1497 xSum *= x; 1498 } 1499 ySum *= y; 1500 } 1501 zSum *= z; 1502 } 1503 tSum *= t; 1504 } 1505 return (sums); 1506 } 1507 # endif /* BuildSums4D */ 1443 1508 1444 1509 /****************************************************************************** … … 1447 1512 in the output vector. This routine works on single-precision polynomials with 1448 1513 double precision data. 1514 1515 XXX EAM : this function is now deprecated: psPolynomial2DEvalVector handles F32 and F64 1449 1516 *****************************************************************************/ 1450 1517 psVector *Polynomial2DEvalVectorD( … … 1658 1725 // xSums look like: 1, x, x^2, ... x^(2n+1) 1659 1726 // Build the B and A data structs. 1727 // XXX EAM : use temp pointers eg vB = B->data.F64 to save redirects 1728 // XXX EAM : this function is only valid for data vectors of F64 1660 1729 for (int k = 0; k < f->n; k++) { 1661 if ((mask != NULL) && mask->data.U8[k])1730 if ((mask != NULL) && (mask->data.U8[k] && maskValue)) { 1662 1731 continue; 1732 } 1663 1733 if (x != NULL) { 1664 1734 xSums = BuildSums1D(xSums, x->data.F64[k], nTerm); … … 1670 1740 wt = 1.0; 1671 1741 } else { 1672 // this should filterfErr == 0 values1673 wt = 1.0 / PS_SQR(fErr->data.F64[k]);1742 // this filters fErr == 0 values 1743 wt = (fErr->data.F64[k] == 0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]); 1674 1744 } 1675 1745 for (int i = 0; i < nTerm; i++) { … … 1686 1756 } 1687 1757 1688 // GaussJordan version 1689 if (0) { 1690 // does the solution in place 1691 psGaussJordan (A, B); 1692 1693 // the first nTerm entries in B correspond directly to the desired 1694 // polynomial coefficients. this is only true for the 1D case 1695 for (int k = 0; k < nTerm; k++) { 1696 myPoly->coeff[k] = B->data.F64[k]; 1697 } 1698 } else { 1699 // LUD version of the fit 1700 psImage *ALUD = NULL; 1701 psVector* outPerm = NULL; 1702 psVector* coeffs = NULL; 1703 1704 ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); 1705 ALUD = psMatrixLUD(ALUD, &outPerm, A); 1706 coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); 1707 for (int k = 0; k < nTerm; k++) { 1708 myPoly->coeff[k] = coeffs->data.F64[k]; 1709 } 1710 psFree(ALUD); 1711 psFree(coeffs); 1712 psFree(outPerm); 1758 // does the solution in place 1759 psGaussJordan (A, B); 1760 1761 // the first nTerm entries in B correspond directly to the desired 1762 // polynomial coefficients. this is only true for the 1D case 1763 for (int k = 0; k < nTerm; k++) { 1764 myPoly->coeff[k] = B->data.F64[k]; 1713 1765 } 1714 1766 … … 1819 1871 } 1820 1872 1821 1873 // This function accepts F32 and F64 input vectors. 1822 1874 psPolynomial1D *psVectorClipFitPolynomial1D( 1823 1875 psPolynomial1D *poly, … … 1827 1879 const psVector *f, 1828 1880 const psVector *fErr, 1829 const psVector *x )1881 const psVector *xIn) 1830 1882 { 1831 1883 // Internal pointers for possibly NULL vectors. 1832 psVector *x 32= NULL;1884 psVector *x = NULL; 1833 1885 1834 1886 PS_ASSERT_POLY_NON_NULL(poly, NULL); … … 1836 1888 PS_ASSERT_PTR_NON_NULL(stats, NULL); 1837 1889 PS_ASSERT_VECTOR_NON_NULL(f, NULL); 1838 PS_ASSERT_VECTOR_TYPE (f, PS_TYPE_F32, NULL);1890 PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); 1839 1891 if (mask != NULL) { 1840 1892 PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); 1841 1893 PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); 1842 1894 } 1843 if (x != NULL) {1844 PS_ASSERT_VECTORS_SIZE_EQUAL(f, x , NULL);1845 PS_ASSERT_VECTOR_TYPE (x, PS_TYPE_F32, NULL);1895 if (xIn != NULL) { 1896 PS_ASSERT_VECTORS_SIZE_EQUAL(f, xIn, NULL); 1897 PS_ASSERT_VECTOR_TYPE_F32_OR_F64(xIn, NULL); 1846 1898 } 1847 1899 if (fErr != NULL) { 1848 1900 PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); 1849 PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); 1850 } 1851 1852 psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); 1901 PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); 1902 } 1903 1904 // assign sequence vector if xIn is NULL 1905 if (xIn == NULL) { 1906 x = psVectorCreate (NULL, 0, f->n, 1, f->type.type); 1907 } else { 1908 x = (psVector *) xIn; 1909 } 1910 1911 // clipping range defined by min and max and/or clipSigma 1912 float minClipSigma; 1913 float maxClipSigma; 1914 if (isfinite(stats->max)) { 1915 maxClipSigma = fabs(stats->clipSigma); 1916 } else { 1917 maxClipSigma = fabs(stats->max); 1918 } 1919 if (isfinite(stats->min)) { 1920 minClipSigma = fabs(stats->clipSigma); 1921 } else { 1922 minClipSigma = fabs(stats->min); 1923 } 1924 psVector *fit = NULL; 1925 psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64); 1926 1927 // eventual expansion: user supplies one of various stats option pairs, 1928 // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to 1929 // evaluate the clipping sigma 1930 // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used 1931 stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); 1932 1933 for (int N = 0; N < stats->clipIter; N++) { 1934 int Nkeep = 0; 1935 1936 poly = psVectorFitPolynomial1D (poly, mask, maskValue, f, fErr, x); 1937 fit = psPolynomial1DEvalVector (poly, x); 1938 resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit); 1939 1940 stats = psVectorStats (stats, resid, NULL, mask, maskValue); 1941 float minClipValue = -minClipSigma*stats->sampleStdev; 1942 float maxClipValue = +maxClipSigma*stats->sampleStdev; 1943 1944 // set mask if pts are not valid 1945 // we are masking out any point which is out of range 1946 // recovery is not allowed with this scheme 1947 for (int i = 0; i < resid->n; i++) { 1948 if ((mask != NULL) && (mask->data.U8[i] & maskValue)) { 1949 continue; 1950 } 1951 if (resid->data.F64[i] - stats->sampleMedian > maxClipValue) { 1952 if (mask != NULL) { 1953 mask->data.U8[i] |= 0x01; 1954 } 1955 continue; 1956 } 1957 if (resid->data.F64[i] - stats->sampleMedian < minClipValue) { 1958 if (mask != NULL) { 1959 mask->data.U8[i] |= 0x01; 1960 } 1961 continue; 1962 } 1963 Nkeep ++; 1964 } 1965 1966 psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n", 1967 Nkeep, x->n); 1968 1969 psFree (fit); 1970 } 1853 1971 // Free psVectors that were created for NULL arguments. 1854 if (x == NULL) { 1855 psFree(x32); 1856 } 1857 return(NULL); 1972 if (xIn == NULL) { 1973 psFree(x); 1974 } 1975 // Free other local temporary variables 1976 psFree (resid); 1977 1978 return (poly); 1858 1979 } 1859 1980 … … 1929 2050 // Build the B and A data structs. 1930 2051 for (int k = 0; k < x->n; k++) { 1931 if ((mask != NULL) && mask->data.U8[k])2052 if ((mask != NULL) && (mask->data.U8[k] & maskValue)) { 1932 2053 continue; 2054 } 2055 1933 2056 Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm); 1934 2057 … … 1936 2059 wt = 1.0; 1937 2060 } else { 1938 // XXX: this should probably by fErr^2 !! 1939 // this should filter fErr == 0 values 1940 // XXX: Why isn't this fErr^2? 1941 wt = 1.0 / fErr->data.F64[k]; 2061 // this filters fErr == 0 values 2062 wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]); 1942 2063 } 1943 2064 … … 1962 2083 1963 2084 // does the solution in place 1964 // XXX: Check return codes!1965 2085 psGaussJordan (A, B); 1966 2086 1967 // XXX: Check return codes! 1968 // ALUD = psMatrixLUD(ALUD, &outPerm, A); 1969 // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); 1970 2087 // select the appropriate solution entries 1971 2088 for (int n = 0; n < nXterm; n++) { 1972 2089 for (int m = 0; m < nYterm; m++) { … … 1985 2102 1986 2103 2104 // XXX EAM : I have implemented a single function to handle the mask/nomask cases 2105 // this function can be deprecated 1987 2106 psPolynomial2D* RobustFit2D_nomask( 1988 2107 psPolynomial2D* poly, … … 2251 2370 PS_ASSERT_PTR_NON_NULL(stats, NULL); 2252 2371 PS_ASSERT_VECTOR_NON_NULL(f, NULL); 2253 PS_ASSERT_VECTOR_TYPE (f, PS_TYPE_F32, NULL);2372 PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); 2254 2373 if (mask != NULL) { 2255 2374 PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); … … 2258 2377 PS_ASSERT_VECTOR_NON_NULL(x, NULL); 2259 2378 PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); 2260 PS_ASSERT_VECTOR_TYPE (x, PS_TYPE_F32, NULL);2379 PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); 2261 2380 PS_ASSERT_VECTOR_NON_NULL(y, NULL); 2262 2381 PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); 2263 PS_ASSERT_VECTOR_TYPE (y, PS_TYPE_F32, NULL);2382 PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); 2264 2383 if (fErr != NULL) { 2265 2384 PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); 2266 PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); 2267 } 2268 2269 if (mask == NULL) { 2270 // XXX: Change argument order. 2271 poly = RobustFit2D_nomask(poly, x, y, f, fErr); 2385 PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); 2386 } 2387 2388 // clipping range defined by min and max and/or clipSigma 2389 float minClipSigma; 2390 float maxClipSigma; 2391 if (isfinite(stats->max)) { 2392 maxClipSigma = fabs(stats->max); 2272 2393 } else { 2273 // XXX: Use maskValue. 2274 // XXX: Change argument order. 2275 poly = RobustFit2D(poly, mask, x, y, f, fErr); 2276 } 2394 maxClipSigma = fabs(stats->clipSigma); 2395 } 2396 if (isfinite(stats->min)) { 2397 minClipSigma = fabs(stats->min); 2398 } else { 2399 minClipSigma = fabs(stats->clipSigma); 2400 } 2401 psVector *fit = NULL; 2402 psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64); 2403 2404 // eventual expansion: user supplies one of various stats option pairs, 2405 // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to 2406 // evaluate the clipping sigma 2407 // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used 2408 stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); 2409 2410 for (int N = 0; N < stats->clipIter; N++) { 2411 int Nkeep = 0; 2412 2413 poly = psVectorFitPolynomial2D (poly, mask, maskValue, f, fErr, x, y); 2414 fit = psPolynomial2DEvalVector (poly, x, y); 2415 resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit); 2416 2417 stats = psVectorStats (stats, resid, NULL, mask, maskValue); 2418 float minClipValue = -minClipSigma*stats->sampleStdev; 2419 float maxClipValue = +maxClipSigma*stats->sampleStdev; 2420 2421 // set mask if pts are not valid 2422 // we are masking out any point which is out of range 2423 // recovery is not allowed with this scheme 2424 for (int i = 0; i < resid->n; i++) { 2425 if ((mask != NULL) && (mask->data.U8[i] & maskValue)) { 2426 continue; 2427 } 2428 if (resid->data.F64[i] - stats->sampleMedian > maxClipValue) { 2429 if (mask != NULL) { 2430 mask->data.U8[i] |= 0x01; 2431 } 2432 continue; 2433 } 2434 if (resid->data.F64[i] - stats->sampleMedian < minClipValue) { 2435 if (mask != NULL) { 2436 mask->data.U8[i] |= 0x01; 2437 } 2438 continue; 2439 } 2440 Nkeep ++; 2441 } 2442 2443 psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n", 2444 Nkeep, x->n); 2445 2446 psFree (fit); 2447 } 2448 // Free local temporary variables 2449 psFree (resid); 2277 2450 2278 2451 if (poly == NULL) { -
trunk/psLib/src/math/psPolynomial.c
r5576 r5624 7 7 * polynomials. It also contains a Gaussian functions. 8 8 * 9 * @version $Revision: 1.13 2$ $Name: not supported by cvs2svn $10 * @date $Date: 2005-11- 22 21:40:40$9 * @version $Revision: 1.133 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2005-11-30 02:00:09 $ 11 11 * 12 12 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 725 725 } 726 726 727 psPolynomial4D* psPolynomial4DAlloc( 728 unsigned int nX, 729 unsigned int nY, 730 unsigned int nZ, 731 unsigned int nT, 732 psPolynomialType type) 727 psPolynomial4D* psPolynomial4DAlloc( unsigned int nX, 728 unsigned int nY, 729 unsigned int nZ, 730 unsigned int nT, 731 psPolynomialType type) 733 732 { 734 733 PS_ASSERT_INT_NONNEGATIVE(nX, NULL); … … 785 784 } 786 785 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 786 psF64 psPolynomial1DEval(const psPolynomial1D* poly, 810 787 psF64 x) … … 824 801 } 825 802 803 // this function must accept F32 and F64 input x vectors 826 804 psVector *psPolynomial1DEvalVector(const psPolynomial1D *poly, 827 805 const psVector *x) … … 829 807 PS_ASSERT_POLY_NON_NULL(poly, NULL); 830 808 PS_ASSERT_VECTOR_NON_NULL(x, NULL); 831 PS_ASSERT_VECTOR_TYPE (x, PS_TYPE_F64, NULL);809 PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); 832 810 833 811 psVector *tmp; 834 812 835 tmp = psVectorAlloc(x->n, PS_TYPE_F64); 836 for (unsigned int i=0;i<x->n;i++) { 837 tmp->data.F64[i] = psPolynomial1DEval(poly, x->data.F64[i]); 838 } 839 813 switch (x->type.type) { 814 case PS_TYPE_F64: 815 tmp = psVectorAlloc(x->n, PS_TYPE_F64); 816 for (unsigned int i=0;i<x->n;i++) { 817 tmp->data.F64[i] = psPolynomial1DEval(poly, x->data.F64[i]); 818 } 819 break; 820 case PS_TYPE_F32: 821 tmp = psVectorAlloc(x->n, PS_TYPE_F32); 822 for (unsigned int i=0;i<x->n;i++) { 823 tmp->data.F32[i] = psPolynomial1DEval(poly, x->data.F32[i]); 824 } 825 break; 826 default: 827 psError(PS_ERR_UNKNOWN, false, "invalid input data type.\n"); 828 return (NULL); 829 } 840 830 return(tmp); 841 831 } … … 859 849 } 860 850 851 // this function must support input data types of F32 and F64 852 // all input vectors data types must match (all F32 or all F64) 861 853 psVector *psPolynomial2DEvalVector(const psPolynomial2D *poly, 862 854 const psVector *x, … … 866 858 PS_ASSERT_POLY_NON_NULL(poly, NULL); 867 859 PS_ASSERT_VECTOR_NON_NULL(x, NULL); 868 PS_ASSERT_VECTOR_TYPE (x, PS_TYPE_F64, NULL);860 PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); 869 861 PS_ASSERT_VECTOR_NON_NULL(y, NULL); 870 PS_ASSERT_VECTOR_TYPE (y, PS_TYPE_F64, NULL);862 PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); 871 863 872 864 psVector *tmp; … … 878 870 } 879 871 880 // Create output vector to return 881 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 882 883 // Evaluate the polynomial at the specified points 884 for (unsigned int i=0; i<vecLen; i++) { 885 tmp->data.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); 886 } 887 872 switch (x->type.type) { 873 case PS_TYPE_F32: 874 if (y->type.type != x->type.type) { 875 psError(PS_ERR_UNKNOWN, true, "type mismatch in data vectors"); 876 return (NULL); 877 } 878 879 // Create output vector to return 880 tmp = psVectorAlloc(vecLen, PS_TYPE_F32); 881 882 // Evaluate the polynomial at the specified points 883 for (unsigned int i=0; i<vecLen; i++) { 884 tmp->data.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[i]); 885 } 886 break; 887 case PS_TYPE_F64: 888 if (y->type.type != x->type.type) { 889 psError(PS_ERR_UNKNOWN, true, "type mismatch in data vectors"); 890 return (NULL); 891 } 892 893 // Create output vector to return 894 tmp = psVectorAlloc(vecLen, PS_TYPE_F64); 895 896 // Evaluate the polynomial at the specified points 897 for (unsigned int i=0; i<vecLen; i++) { 898 tmp->data.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); 899 } 900 break; 901 default: 902 psError(PS_ERR_UNKNOWN, false, "invalid input data type.\n"); 903 return (NULL); 904 } 888 905 // Return output vector 889 906 return(tmp); -
trunk/psLib/src/math/psStats.c
r5576 r5624 7 7 * on those data structures. 8 8 * 9 * @author GLG, MHPCC10 *11 * XXX: The following stats members are never used, or set in this code.12 * stats->robustN5013 * stats->clippedNvalues14 * stats->binsize15 *16 *17 *18 *19 * @version $Revision: 1.154$ $Name: not supported by cvs2svn $20 * @date $Date: 2005-11-22 21:40:40$21 *22 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii23 */9 * @author GLG, MHPCC 10 * 11 * XXX: The following stats members are never used, or set in this code. 12 * stats->robustN50 13 * stats->clippedNvalues 14 * stats->binsize 15 * 16 * 17 * 18 * 19 * @version $Revision: 1.155 $ $Name: not supported by cvs2svn $ 20 * @date $Date: 2005-11-30 02:00:09 $ 21 * 22 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 23 */ 24 24 25 25 #include <stdlib.h> … … 1131 1131 statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN); 1132 1132 p_psMemSetPersistent(statsTmp, true); 1133 } else { 1134 // XXX EAM : initialize structure if already allocated 1135 statsTmp->sampleMean = NAN; 1136 statsTmp->sampleMedian = NAN; 1137 statsTmp->sampleStdev = NAN; 1138 statsTmp->sampleUQ = NAN; 1139 statsTmp->sampleLQ = NAN; 1140 statsTmp->robustMean = NAN; 1141 statsTmp->robustMedian = NAN; 1142 statsTmp->robustMode = NAN; 1143 statsTmp->robustStdev = NAN; 1144 statsTmp->robustUQ = NAN; 1145 statsTmp->robustLQ = NAN; 1146 statsTmp->robustN50 = -1; // XXX: This is never used 1147 statsTmp->robustNfit = -1; 1148 statsTmp->clippedMean = NAN; 1149 statsTmp->clippedStdev = NAN; 1150 statsTmp->clippedNvalues = -1; // XXX: This is never used 1151 statsTmp->clipSigma = 3.0; 1152 statsTmp->clipIter = 3; 1153 statsTmp->min = NAN; 1154 statsTmp->max = NAN; 1155 statsTmp->binsize = NAN; // XXX: This is never used 1133 1156 } 1134 1157 -
trunk/psLib/src/mathtypes/psImage.c
r5572 r5624 9 9 * @author Ross Harman, MHPCC 10 10 * 11 * @version $Revision: 1.9 0$ $Name: not supported by cvs2svn $12 * @date $Date: 2005-11- 22 19:58:16$11 * @version $Revision: 1.91 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2005-11-30 02:00:10 $ 13 13 * 14 14 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -
trunk/psLib/test/sys/verified/tst_psMemory.stderr
r5216 r5624 169 169 <HOST>|E|memCheckTypes (FILE:LINENO) 170 170 psMemCheckBitSet failed in memCheckType. 171 <HOST>|E|psPolynomial4DAlloc (FILE:LINENO)172 Error: nX is 0 or less.173 <HOST>|E|psPolynomial4DAlloc (FILE:LINENO)174 Error: nX is 0 or less.175 <HOST>|E|psPolynomial2DAlloc (FILE:LINENO)176 Error: nX is 0 or less.177 <HOST>|E|psPolynomial2DAlloc (FILE:LINENO)178 Error: nX is 0 or less.179 171 180 172 ---> TESTPOINT PASSED (psMemory{psMemCheckType} | tst_psMemory.c)
Note:
See TracChangeset
for help on using the changeset viewer.
