Changeset 3514
- Timestamp:
- Mar 25, 2005, 5:44:31 PM (21 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/pmObjects.c (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/pmObjects.c
r3512 r3514 5 5 * @author GLG, MHPCC 6 6 * 7 * @version $Revision: 1. 9$ $Name: not supported by cvs2svn $8 * @date $Date: 2005-03-2 5 23:23:29$7 * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2005-03-26 03:44:31 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 1082 1082 XXX: This should probably be a public function. 1083 1083 1084 XXX: We don't care about the model derivatives here. So, for the individual1085 model evaluation functions, we should probably separate the derivative1086 calculation.1087 1088 XXX: Get rid of deriv.1089 1090 1084 XXX: Use static vectors for x. 1091 1085 … … 1096 1090 psU32 col) 1097 1091 { 1098 // XXX: You won't need this once you separate the function evaluation.1099 // from the derivative calculation.1100 psVector *deriv = psVectorAlloc(src->models->Nparams, PS_TYPE_F32);1101 1092 // XXX: The following step will not be necessary if the models->params 1102 1093 // member is a psVector. … … 1113 1104 switch (src->models->type) { 1114 1105 case PS_MODEL_GAUSS: 1115 tmpF = pmMinLM_Gauss2D( deriv, params, x);1106 tmpF = pmMinLM_Gauss2D(NULL, params, x); 1116 1107 break; 1117 1108 case PS_MODEL_PGAUSS: 1118 tmpF = pmMinLM_PsuedoGauss2D( deriv, params, x);1109 tmpF = pmMinLM_PsuedoGauss2D(NULL, params, x); 1119 1110 break; 1120 1111 case PS_MODEL_TWIST_GAUSS: 1121 tmpF = pmMinLM_TwistGauss2D( deriv, params, x);1112 tmpF = pmMinLM_TwistGauss2D(NULL, params, x); 1122 1113 break; 1123 1114 case PS_MODEL_WAUSS: 1124 tmpF = pmMinLM_Wauss2D( deriv, params, x);1115 tmpF = pmMinLM_Wauss2D(NULL, params, x); 1125 1116 break; 1126 1117 case PS_MODEL_SERSIC: 1127 tmpF = pmMinLM_Sersic( deriv, params, x);1118 tmpF = pmMinLM_Sersic(NULL, params, x); 1128 1119 break; 1129 1120 case PS_MODEL_SERSIC_CORE: 1130 tmpF = pmMinLM_SersicCore( deriv, params, x);1121 tmpF = pmMinLM_SersicCore(NULL, params, x); 1131 1122 break; 1132 1123 default: … … 1135 1126 } 1136 1127 1137 psFree(deriv);1138 1128 psFree(params); 1139 1129 psFree(x); … … 1388 1378 #define MAX_PARAMS 100 1389 1379 1390 psVector *deriv = psVectorAlloc(MAX_PARAMS, PS_TYPE_F32);1391 1380 psVector *params = psVectorAlloc(MAX_PARAMS, PS_TYPE_F32); 1392 1381 psVector *x = psVectorAlloc(2, PS_TYPE_F32); … … 1410 1399 switch (src->models->type) { 1411 1400 case PS_MODEL_GAUSS: 1412 pixelValue = pmMinLM_Gauss2D( deriv, params, x);1401 pixelValue = pmMinLM_Gauss2D(NULL, params, x); 1413 1402 break; 1414 1403 case PS_MODEL_PGAUSS: 1415 pixelValue = pmMinLM_PsuedoGauss2D( deriv, params, x);1404 pixelValue = pmMinLM_PsuedoGauss2D(NULL, params, x); 1416 1405 break; 1417 1406 case PS_MODEL_TWIST_GAUSS: 1418 pixelValue = pmMinLM_TwistGauss2D( deriv, params, x);1407 pixelValue = pmMinLM_TwistGauss2D(NULL, params, x); 1419 1408 break; 1420 1409 case PS_MODEL_WAUSS: 1421 pixelValue = pmMinLM_Wauss2D( deriv, params, x);1410 pixelValue = pmMinLM_Wauss2D(NULL, params, x); 1422 1411 break; 1423 1412 case PS_MODEL_SERSIC: 1424 pixelValue = pmMinLM_Sersic( deriv, params, x);1413 pixelValue = pmMinLM_Sersic(NULL, params, x); 1425 1414 break; 1426 1415 case PS_MODEL_SERSIC_CORE: 1427 pixelValue = pmMinLM_SersicCore( deriv, params, x);1416 pixelValue = pmMinLM_SersicCore(NULL, params, x); 1428 1417 break; 1429 1418 default: 1430 1419 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1431 1420 psFree(x); 1432 psFree(deriv);1433 1421 psFree(params); 1434 1422 return(false); … … 1447 1435 } 1448 1436 psFree(x); 1449 psFree(deriv);1450 1437 psFree(params); 1451 1438 return(true); … … 1499 1486 XXX: Consider getting rid of the parameter checks since this might consume 1500 1487 a significant fraction of this function CPU time. 1488 1489 XXX: I added the following. Must conform with IfA. If deriv==NULL, then 1490 we simply don't perform the derivative calculations. 1501 1491 *****************************************************************************/ 1502 1492 psF32 pmMinLM_Gauss2D(psVector *deriv, … … 1504 1494 psVector *x) 1505 1495 { 1506 PS_VECTOR_CHECK_NULL(deriv, NAN);1507 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);1508 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);1509 1496 PS_VECTOR_CHECK_NULL(params, NAN); 1510 1497 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); … … 1522 1509 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1523 1510 1524 psF32 q = params->data.F32[1]*r; 1525 deriv->data.F32[0] = +1.0; 1526 deriv->data.F32[1] = +r; 1527 deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); 1528 deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); 1529 deriv->data.F32[4] = -2.0*q*px*X; 1530 deriv->data.F32[5] = -2.0*q*py*Y; 1531 deriv->data.F32[6] = -q*X*Y; 1511 if (deriv != NULL) { 1512 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1513 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); 1514 1515 psF32 q = params->data.F32[1]*r; 1516 deriv->data.F32[0] = +1.0; 1517 deriv->data.F32[1] = +r; 1518 deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); 1519 deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); 1520 deriv->data.F32[4] = -2.0*q*px*X; 1521 deriv->data.F32[5] = -2.0*q*py*Y; 1522 deriv->data.F32[6] = -q*X*Y; 1523 } 1532 1524 1533 1525 return(f); … … 1588 1580 psVector *x) 1589 1581 { 1590 PS_VECTOR_CHECK_NULL(deriv, NAN);1591 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);1592 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);1593 1582 PS_VECTOR_CHECK_NULL(params, NAN); 1594 1583 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); … … 1607 1596 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1608 1597 1609 // 1610 // note difference from a pure gaussian: q = params->data.F32[1]*r 1611 // 1612 1613 psF32 q = params->data.F32[1]*r*r*t; 1614 deriv->data.F32[0] = +1.0; 1615 deriv->data.F32[1] = +r; 1616 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1617 deriv->data.F32[3] = q * 1618 (2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1619 deriv->data.F32[4] = -2.0*q*px*X; 1620 deriv->data.F32[5] = -2.0*q*py*Y; 1621 deriv->data.F32[6] = -q*X*Y; 1598 if (deriv != NULL) { 1599 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1600 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); 1601 1602 // 1603 // note difference from a pure gaussian: q = params->data.F32[1]*r 1604 // 1605 psF32 q = params->data.F32[1]*r*r*t; 1606 deriv->data.F32[0] = +1.0; 1607 deriv->data.F32[1] = +r; 1608 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1609 deriv->data.F32[3] = q * 1610 (2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1611 deriv->data.F32[4] = -2.0*q*px*X; 1612 deriv->data.F32[5] = -2.0*q*py*Y; 1613 deriv->data.F32[6] = -q*X*Y; 1614 } 1622 1615 1623 1616 return(f); … … 1682 1675 psVector *x) 1683 1676 { 1684 PS_VECTOR_CHECK_NULL(deriv, NAN);1685 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);1686 PS_VECTOR_CHECK_SIZE(deriv, 9, NAN);1687 1677 PS_VECTOR_CHECK_NULL(params, NAN); 1688 1678 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); … … 1701 1691 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1702 1692 1703 // 1704 // note difference from gaussian: q = params->data.F32[1]*r 1705 // 1706 psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); 1707 deriv->data.F32[0] = +1.0; 1708 deriv->data.F32[1] = +r; 1709 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1710 deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1711 deriv->data.F32[4] = -2.0*q*px*X; 1712 deriv->data.F32[5] = -2.0*q*py*Y; 1713 deriv->data.F32[6] = -q*X*Y; 1714 deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; 1715 deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; 1716 // 1717 // The values of 100 dampen the swing of params->data.F32[7,8] */ 1718 // 1693 if (deriv != NULL) { 1694 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1695 PS_VECTOR_CHECK_SIZE(deriv, 9, NAN); 1696 // 1697 // note difference from gaussian: q = params->data.F32[1]*r 1698 // 1699 1700 psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); 1701 deriv->data.F32[0] = +1.0; 1702 deriv->data.F32[1] = +r; 1703 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1704 deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1705 deriv->data.F32[4] = -2.0*q*px*X; 1706 deriv->data.F32[5] = -2.0*q*py*Y; 1707 deriv->data.F32[6] = -q*X*Y; 1708 deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; 1709 deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; 1710 // 1711 // The values of 100 dampen the swing of params->data.F32[7,8] */ 1712 // 1713 } 1719 1714 1720 1715 return(f); … … 1786 1781 psVector *x) 1787 1782 { 1788 PS_VECTOR_CHECK_NULL(deriv, NAN);1789 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);1790 PS_VECTOR_CHECK_SIZE(deriv, 11, NAN);1791 1783 PS_VECTOR_CHECK_NULL(params, NAN); 1792 1784 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); … … 1808 1800 1809 1801 psF32 f = params->data.F32[5]*r + params->data.F32[6]; 1810 psF32 q1 = params->data.F32[5]*PS_SQR(r); 1811 psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); 1812 deriv->data.F32[0] = +1.0; 1813 deriv->data.F32[1] = +r; 1814 deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y); 1815 deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X); 1816 1817 // 1818 // These fudge factors impede the growth of params->data.F32[4] beyond 1819 // params->data.F32[7]. 1820 // 1821 psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); 1822 psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 1823 deriv->data.F32[4] = -2.0*q1*px1*X*f2; 1824 1825 // 1826 // These fudge factors impede the growth of params->data.F32[5] beyond 1827 // params->data.F32[8]. 1828 // 1829 f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); 1830 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 1831 deriv->data.F32[5] = -2.0*q1*py1*Y*f2; 1832 deriv->data.F32[6] = -q1*X*Y; 1833 deriv->data.F32[7] = -2.0*q2*px2*X; 1834 deriv->data.F32[8] = -2.0*q2*py2*Y; 1835 deriv->data.F32[9] = -q2*X*Y; 1836 deriv->data.F32[10] = -q1*log(z2); 1802 1803 if (deriv != NULL) { 1804 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1805 PS_VECTOR_CHECK_SIZE(deriv, 11, NAN); 1806 1807 psF32 q1 = params->data.F32[5]*PS_SQR(r); 1808 psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); 1809 deriv->data.F32[0] = +1.0; 1810 deriv->data.F32[1] = +r; 1811 deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y); 1812 deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X); 1813 1814 // 1815 // These fudge factors impede the growth of params->data.F32[4] beyond 1816 // params->data.F32[7]. 1817 // 1818 psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); 1819 psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 1820 deriv->data.F32[4] = -2.0*q1*px1*X*f2; 1821 1822 // 1823 // These fudge factors impede the growth of params->data.F32[5] beyond 1824 // params->data.F32[8]. 1825 // 1826 f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); 1827 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 1828 deriv->data.F32[5] = -2.0*q1*py1*Y*f2; 1829 deriv->data.F32[6] = -q1*X*Y; 1830 deriv->data.F32[7] = -2.0*q2*px2*X; 1831 deriv->data.F32[8] = -2.0*q2*py2*Y; 1832 deriv->data.F32[9] = -q2*X*Y; 1833 deriv->data.F32[10] = -q1*log(z2); 1834 } 1837 1835 1838 1836 return(f); … … 1896 1894 psVector *x) 1897 1895 { 1898 PS_VECTOR_CHECK_NULL(deriv, NAN);1899 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);1900 PS_VECTOR_CHECK_SIZE(deriv, 8, NAN);1901 1896 PS_VECTOR_CHECK_NULL(params, NAN); 1902 1897 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); … … 1905 1900 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1906 1901 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1902 1903 if (deriv != NULL) { 1904 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1905 PS_VECTOR_CHECK_SIZE(deriv, 8, NAN); 1906 } 1907 1907 1908 1908 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); … … 1968 1968 psVector *x) 1969 1969 { 1970 PS_VECTOR_CHECK_NULL(deriv, NAN);1971 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);1972 PS_VECTOR_CHECK_SIZE(deriv, 12, NAN);1973 1970 PS_VECTOR_CHECK_NULL(params, NAN); 1974 1971 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); … … 1977 1974 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1978 1975 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1976 1977 if (deriv != NULL) { 1978 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1979 PS_VECTOR_CHECK_SIZE(deriv, 12, NAN); 1980 } 1979 1981 1980 1982 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
Note:
See TracChangeset
for help on using the changeset viewer.
