Index: /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c (revision 36318) +++ /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c (revision 36319) @@ -27,14 +27,29 @@ // fixed set of radii and number of entries in circularly symmetric profile -# define NRAD_MAX 44 -static float radii2[NRAD_MAX] = { 0.0, 1.0, 2.0, 4.0, 5.0, - 8.0, 9.0, 10.0, 13.0, 16.0, - 17.0, 18.0, 20.0, 25.0, 26.0, - 29.0, 32.0, 34.0, 36.0, 37.0, - 40.0, 41.0, 45.0, 49.0, 50.0, - 52.0, 53.0, 58.0, 61.0, 64.0, - 65.0, 68.0, 72.0, 73.0, 74.0, - 80.0, 81.0, 82.0, 85.0, 89.0, - 90.0, 97.0, 98.0, 100.0 +# define NRAD_MAX 87 +# define RAD_MAX 15 +static float radii2[NRAD_MAX] = { 0.0, 1.0, 2.0, 4.0, 5.0, + 8.0, 9.0, 10.0, 13.0, 16.0, + 17.0, 18.0, 20.0, 25.0, 26.0, + + 29.0, 32.0, 34.0, 36.0, 37.0, + 40.0, 41.0, 45.0, 49.0, 50.0, + 52.0, 53.0, 58.0, 61.0, 64.0, + + 65.0, 68.0, 72.0, 73.0, 74.0, + 80.0, 81.0, 82.0, 85.0, 89.0, + 90.0, 97.0, 98.0, 100.0, 101.0, + + 104.0, 106.0, 109.0, 113.0, 116.0, + 117.0, 121.0, 122.0, 125.0, 128.0, + 130.0, 136.0, 137.0, 144.0, 145.0, + + 146.0, 148.0, 149.0, 153.0, 157.0, + 160.0, 162.0, 164.0, 169.0, 170.0, + 173.0, 178.0, 180.0, 181.0, 185.0, + + 193.0, 194.0, 196.0, 197.0, 200.0, + 202.0, 205.0, 208.0, 212.0, 218.0, + 221.0, 225.0 }; @@ -42,10 +57,24 @@ 4, 4, 8, 8, 4, 8, 4, 8, 12, 8, + 8, 4, 8, 4, 8, 8, 8, 8, 4, 12, 8, 8, 8, 8, 4, + 16, 8, 4, 8, 8, - 8, 4, 8, 16, 8, - 8, 8, 4, 12 + 8, 4, 8, 16, 8, + 8, 8, 4, 12, 8, + + 8, 8, 8, 8, 8, + 8, 4, 8, 16, 4, + 16, 8, 8, 8, 16, + + 8, 8, 8, 8, 8, + 8, 4, 8, 12, 16, + 8, 8, 8, 8, 16, + + 8, 8, 4, 8, 12, + 8, 16, 8, 8, 8, + 16, 12 }; @@ -85,4 +114,6 @@ psImageSmooth2dCacheData *psImageSmooth2dCacheAlloc (float Nsigma) { + psAssert (isfinite(Nsigma) && (Nsigma > 0.5), "Nsigma is not valid"); + psImageSmooth2dCacheData *smdata = psAlloc(sizeof(psImageSmooth2dCacheData)); psMemSetDeallocator(smdata, (psFreeFunc) psImageSmooth2dCacheDataFree); @@ -100,5 +131,5 @@ // check for NULL structure elements? int Ns = (int)(smdata->Nsigma * sigma); - Ns = PS_MAX (3, PS_MIN (Ns, 10)); + Ns = PS_MAX (3, PS_MIN (Ns, RAD_MAX)); smdata->Ns = Ns; @@ -108,5 +139,5 @@ smdata->radflux = psAlloc(sizeof(float)*NRAD_MAX); - float sum = 0.0; + double sum = 0.0; for (int i = 0; i < NRAD_MAX; i++) { if (radii2[i] > Ns2) { @@ -130,5 +161,5 @@ // check for NULL structure elements? int Ns = (int)(smdata->Nsigma * sigma); - Ns = PS_MAX (3, PS_MIN (Ns, 10)); + Ns = PS_MAX (3, PS_MIN (Ns, RAD_MAX)); smdata->Ns = Ns; @@ -275,4 +306,67 @@ if (Ns <= 10) goto finish; + ADD_RAND ( 44, 1, 10); // r^2 = 101 + ADD_RAND ( 45, 2, 10); // r^2 = 104 + ADD_RAND ( 46, 5, 9); // r^2 = 106 + ADD_RAND ( 47, 3, 10); // r^2 = 109 + ADD_RAND ( 48, 7, 8); // r^2 = 113 + ADD_RAND ( 49, 4, 10); // r^2 = 116 + ADD_RAND ( 50, 6, 9); // r^2 = 117 + ADD_AXIS ( 51, 11) ; // r^2 = 121 + if (Ns <= 11) goto finish; + + ADD_RAND ( 52, 1, 11); // r^2 = 122 + ADD_RAND ( 53, 2, 11); // r^2 = 125 * + ADD_RAND ( 53, 5, 10); // r^2 = 125 * + ADD_DIAG ( 54, 8) ; // r^2 = 128 + ADD_RAND ( 55, 3, 11); // r^2 = 130 * + ADD_RAND ( 55, 7, 9); // r^2 = 130 * + ADD_RAND ( 56, 6, 10); // r^2 = 136 + ADD_RAND ( 57, 4, 11); // r^2 = 137 + ADD_AXIS ( 58, 12) ; // r^2 = 144 + if (Ns <= 12) goto finish; + + ADD_RAND ( 59, 1, 12); // r^2 = 145 * + ADD_RAND ( 59, 8, 9); // r^2 = 145 * + ADD_RAND ( 60, 5, 11); // r^2 = 146 + ADD_RAND ( 61, 2, 12); // r^2 = 148 + ADD_RAND ( 62, 7, 10); // r^2 = 149 + ADD_RAND ( 63, 3, 12); // r^2 = 153 + ADD_RAND ( 64, 6, 11); // r^2 = 157 + ADD_RAND ( 65, 4, 12); // r^2 = 160 + ADD_DIAG ( 66, 9) ; // r^2 = 162 + ADD_RAND ( 67, 8, 10); // r^2 = 164 + ADD_RAND ( 68, 5, 12); // r^2 = 169 * + ADD_AXIS ( 68, 13) ; // r^2 = 169 * + if (Ns <= 13) goto finish; + + ADD_RAND ( 69, 1, 13); // r^2 = 170 * + ADD_RAND ( 69, 7, 11); // r^2 = 170 * + ADD_RAND ( 70, 2, 13); // r^2 = 173 + ADD_RAND ( 71, 3, 13); // r^2 = 178 + ADD_RAND ( 72, 6, 12); // r^2 = 180 + ADD_RAND ( 73, 9, 10); // r^2 = 181 + ADD_RAND ( 74, 4, 13); // r^2 = 185 * + ADD_RAND ( 74, 8, 11); // r^2 = 185 * + ADD_RAND ( 75, 7, 12); // r^2 = 193 + ADD_RAND ( 76, 5, 13); // r^2 = 194 + ADD_AXIS ( 77, 14) ; // r^2 = 196 + if (Ns <= 14) goto finish; + + ADD_RAND ( 78, 1, 14); // r^2 = 197 + ADD_DIAG ( 79, 10) ; // r^2 = 200 * + ADD_RAND ( 79, 2, 14); // r^2 = 200 * + ADD_RAND ( 80, 9, 11); // r^2 = 202 + ADD_RAND ( 81, 3, 14); // r^2 = 205 * + ADD_RAND ( 81, 6, 13); // r^2 = 205 * + ADD_RAND ( 82, 8, 12); // r^2 = 208 + ADD_RAND ( 83, 4, 14); // r^2 = 212 + ADD_RAND ( 84, 7, 13); // r^2 = 218 + ADD_RAND ( 85, 10, 11); // r^2 = 221 * + ADD_RAND ( 85, 5, 14); // r^2 = 221 * + ADD_RAND ( 86, 9, 12); // r^2 = 225 * + ADD_AXIS ( 86, 15); // r^2 = 225 * + if (Ns <= 15) goto finish; + finish: vo[iy][ix] = s; @@ -298,2 +392,14 @@ } +/* +101, 104, 106, 109, 113, +116, 117, 121, 122, 125, +128, 130, 136, 137, 144, +145, 146, 148, 149, 153, +157, 160, 162, 164, 169, +170, 173, 178, 180, 181, +185, 193, 194, 196, 197, +200, 202, 205, 208, 212, +218, 221, 225 +*/ + Index: /branches/eam_branches/ipp-20130904/psLib/test/imageops/tap_psImageConvolve2dCache.c =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/test/imageops/tap_psImageConvolve2dCache.c (revision 36318) +++ /branches/eam_branches/ipp-20130904/psLib/test/imageops/tap_psImageConvolve2dCache.c (revision 36319) @@ -14,5 +14,5 @@ #include "pstap.h" -bool tap_convolve_delta (int Nx, int Ny, int Xo, int Yo, float sigma, float kappa); +bool tap_convolve_delta (int Nx, int Ny, int Xo, int Yo, float Nsigma, float sigma, float kappa); int main(int argc, char **argv) { @@ -23,27 +23,45 @@ // generate an image to convolve - tap_convolve_delta (301, 301, 150, 150, 2.0, 1.0); - tap_convolve_delta (301, 301, 100, 200, 2.0, 1.0); - tap_convolve_delta (301, 301, 200, 100, 2.0, 1.0); - tap_convolve_delta (301, 301, 200, 200, 2.0, 1.0); - tap_convolve_delta (301, 301, 100, 100, 2.0, 1.0); + tap_convolve_delta (301, 301, 150, 150, 3.0, 2.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 3.0, 2.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 3.0, 2.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 3.0, 2.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 3.0, 2.0, 1.0); - tap_convolve_delta (301, 301, 150, 150, 3.0, 1.0); - tap_convolve_delta (301, 301, 100, 200, 3.0, 1.0); - tap_convolve_delta (301, 301, 200, 100, 3.0, 1.0); - tap_convolve_delta (301, 301, 200, 200, 3.0, 1.0); - tap_convolve_delta (301, 301, 100, 100, 3.0, 1.0); + tap_convolve_delta (301, 301, 150, 150, 3.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 3.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 3.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 3.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 3.0, 3.0, 1.0); - tap_convolve_delta (301, 301, 150, 150, 1.0, 1.0); - tap_convolve_delta (301, 301, 100, 200, 1.0, 1.0); - tap_convolve_delta (301, 301, 200, 100, 1.0, 1.0); - tap_convolve_delta (301, 301, 200, 200, 1.0, 1.0); - tap_convolve_delta (301, 301, 100, 100, 1.0, 1.0); + tap_convolve_delta (301, 301, 150, 150, 3.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 3.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 3.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 3.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 3.0, 1.0, 1.0); - tap_convolve_delta (301, 301, 150, 150, 2.0, 2.0); - tap_convolve_delta (301, 301, 100, 200, 2.0, 2.0); - tap_convolve_delta (301, 301, 200, 100, 2.0, 2.0); - tap_convolve_delta (301, 301, 200, 200, 2.0, 2.0); - tap_convolve_delta (301, 301, 100, 100, 2.0, 2.0); + tap_convolve_delta (301, 301, 150, 150, 3.0, 2.0, 2.0); + tap_convolve_delta (301, 301, 100, 200, 3.0, 2.0, 2.0); + tap_convolve_delta (301, 301, 200, 100, 3.0, 2.0, 2.0); + tap_convolve_delta (301, 301, 200, 200, 3.0, 2.0, 2.0); + tap_convolve_delta (301, 301, 100, 100, 3.0, 2.0, 2.0); + + tap_convolve_delta (301, 301, 150, 150, 4.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 4.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 4.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 4.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 4.0, 3.0, 1.0); + + tap_convolve_delta (301, 301, 150, 150, 5.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 5.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 5.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 5.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 5.0, 3.0, 1.0); + + tap_convolve_delta (301, 301, 150, 150, 5.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 5.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 5.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 5.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 5.0, 1.0, 1.0); ok(!psMemCheckLeaks (id, NULL, stdout, false), "no memory leaks"); @@ -52,5 +70,5 @@ } -bool tap_convolve_delta (int Nx, int Ny, int Xo, int Yo, float sigma, float kappa) { +bool tap_convolve_delta (int Nx, int Ny, int Xo, int Yo, float Nsigma, float sigma, float kappa) { bool status; @@ -60,8 +78,8 @@ input->data.F32[Yo][Xo] = 1.0; - psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (3.0); + psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (Nsigma); ok(smdata, "init smdata"); ok(!smdata->radflux, "init radflux to NULL"); - ok(smdata->Nsigma == 3, "init Nsigma to 3"); + ok(smdata->Nsigma == Nsigma, "init Nsigma"); status = psImageSmooth2dCacheKernel_PS1_V1 (smdata, sigma, kappa); @@ -71,5 +89,5 @@ ok(status, "smoothed image"); - float s = 0.0; + double s = 0.0; psImage *model = psImageAlloc (Nx, Ny, PS_TYPE_F32); @@ -77,6 +95,6 @@ for (int iy = 0; iy < model->numRows; iy++) { for (int ix = 0; ix < model->numCols; ix++) { - float z = PS_SQR((ix + 0.5 - Xo - 0.5) / sigma) + PS_SQR((iy + 0.5 - Yo - 0.5) / sigma); - if (z > PS_SQR(smdata->Nsigma)) continue; + float z = 0.5*PS_SQR((ix + 0.5 - Xo - 0.5) / sigma) + 0.5*PS_SQR((iy + 0.5 - Yo - 0.5) / sigma); + if (z > 0.5*PS_SQR(smdata->Nsigma)) continue; float f = 1.0 / (1.0 + z*kappa + pow(z, 1.666)); @@ -101,8 +119,8 @@ } } - ok (fabs(min_df) < 1e-5, "minor deviations"); - ok (fabs(max_df) < 1e-5, "minor deviations"); + ok (fabs(min_df) < 3e-5, "minor deviations"); + ok (fabs(max_df) < 3e-5, "minor deviations"); - if ((fabs(min_df) >= 1e-5) || (fabs(max_df) >= 1e-5)) { + if ((fabs(min_df) >= 3e-5) || (fabs(max_df) >= 3e-5)) { psFits *fits = NULL; fits = psFitsOpen ("input.fits", "w");