From 46dfd1884aee0faffa3d5f95f56bc3deab4823d9 Mon Sep 17 00:00:00 2001 From: GillesDuvert Date: Sat, 30 Nov 2024 19:24:13 +0100 Subject: [PATCH 1/2] solved #1906 --- src/widget.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/widget.cpp b/src/widget.cpp index 94325fed9..336bba135 100644 --- a/src/widget.cpp +++ b/src/widget.cpp @@ -2511,9 +2511,7 @@ BaseGDL* widget_info( EnvT* e ) { id = (*static_cast (ev->GetTag(idIx, 0)))[0]; // get its id for (SizeT i = 0; i < widgetIDList.size(); i++) { //is ID corresponding to any widget in list? if (widgetIDList.at(i) == id) { //if yes - //IMPORTANT: return ev immediately. This is what permits #1685: an event trapped by WIDGET_EVENT does not behave - //like the same event processed in the evenloop. - return ev; + goto endwait; } } GDLWidget::widgetEventQueue.PushBack(ev); From 2dce597b380f050aa8f1b72afc64bd9bb37b1728 Mon Sep 17 00:00:00 2001 From: GillesDuvert Date: Tue, 3 Dec 2024 12:16:57 +0100 Subject: [PATCH 2/2] added support for COMPLEX values in POLY_2D slight rewriting of code for simplification. Simple tests pass. --- src/poly_2d.cpp | 523 +++++++++++------- ...poly2d_complex_handle_cubic_function.incpp | 18 + .../poly2d_complex_handle_function.incpp | 18 + testsuite/test_all_basic_functions.pro | 2 +- testsuite/test_bugs_poly2d.pro | 37 +- 5 files changed, 383 insertions(+), 215 deletions(-) create mode 100644 src/snippets/poly2d_complex_handle_cubic_function.incpp create mode 100644 src/snippets/poly2d_complex_handle_function.incpp diff --git a/src/poly_2d.cpp b/src/poly_2d.cpp index 391fbc6fa..eed115961 100644 --- a/src/poly_2d.cpp +++ b/src/poly_2d.cpp @@ -76,16 +76,12 @@ namespace lib { } - DFloat * poly2d_compute_init_x(poly2d * p, SizeT n) { + DFloat * poly2d_compute_init(poly2d * p, const SizeT m0, const SizeT m1) { + SizeT n=(m1>m0)?m1:m0; DFloat * res= (DFloat*) malloc(p->nc*n*sizeof(DFloat)); for (auto i=0, s=0; i< n; ++i) for (DLong k = 0; k < p->nc; k++) res[s++]=ipowI(i, p->px[k]); return res; } - DFloat* poly2d_compute_init_y(poly2d * p, SizeT n) { - DFloat * res= (DFloat*) malloc(p->nc*n*sizeof(DFloat)); - for (auto i=0,s=0; i< n; ++i) for (DLong k = 0; k < p->nc; k++) res[s++]=ipowI(i, p->py[k]); - return res; - } // cubic kernels are much faster using a precoputed table of kernel values digitzed on 1/1000 th of a pixel. static DFloat* cubicKernel; @@ -126,6 +122,9 @@ DFloat bicubicInterpolate (DFloat p[4][4], DFloat x, DFloat y, DFloat * precompu The returned array of DFloats must be deallocated using free(). */ /*--------------------------------------------------------------------------*/ + +// commenting unused code is the best way to avoid a frustrating uncompleteness of coverage +// when using coverage tests. DFloat * generate_interpolation_kernel(/* int kernel_type, */ DFloat cubicParameter) { DFloat * tab ; @@ -196,8 +195,8 @@ template< typename T1, typename T2> const SizeT nCols, const SizeT nRows, BaseGDL* data_, - DFloat * const P, - DFloat * const Q, + const DFloat * P, + const DFloat * Q, DFloat const initvalue_, const bool doMissing) { // std::cerr<<"warp_linear0\n"; @@ -335,14 +334,14 @@ template< typename T1, typename T2> BaseGDL* warp0( const SizeT nCols, const SizeT nRows, - BaseGDL* const data_, - poly2d* const poly_u, - poly2d* const poly_v, + BaseGDL * data_, + poly2d * poly_u, + poly2d * poly_v, const DFloat initvalue_, - bool doMissing) { + const bool doMissing) { // std::cerr<<"warp0\n"; - const SizeT lx = data_->Dim(0); - const SizeT ly = data_->Dim(1); + SizeT lx = data_->Dim(0); + SizeT ly = data_->Dim(1); dimension dim(nCols, nRows); T1* res_ = new T1(dim, BaseGDL::NOZERO); @@ -366,19 +365,17 @@ template< typename T1, typename T2> DFloat fl0y=0; //these are accelerators - not using them increase exec time by 2 or 3. DLong nc=poly_u->nc; - DFloat * const xcoefu=poly2d_compute_init_x(poly_u,lx); - DFloat * const ycoefu=poly2d_compute_init_y(poly_u,lx); - DFloat * const xcoefv=poly2d_compute_init_x(poly_v,ly); - DFloat * const ycoefv=poly2d_compute_init_y(poly_v,ly); + assert(poly_u->nc == poly_v->nc); + DFloat * const coef=poly2d_compute_init(poly_u,lx,ly); //poly_u and poly_va share the same size. /* Double loop on the output image */ if (doMissing) { if ((GDL_NTHREADS = parallelize(nEl)) == 1) { for (OMPInt j = 0; j < nRows; ++j) { for (OMPInt i = 0; i < nCols; ++i) { - DFloat x = poly_u->c[0]; for (auto k=1; k< nc; ++k) x+=poly_u->c[k]*xcoefu[j*nc+k]*ycoefu[i*nc+k]; + DFloat x = poly_u->c[0]; for (auto k=1; k< nc; ++k) x+=poly_u->c[k]*coef[j*nc+k]*coef[i*nc+k]; if (x < fl0x) continue; if (x >= fllx) continue; // already initialised to 'missing' value. - DFloat y = poly_v->c[0]; for (auto k=1; k< nc; ++k) y+=poly_v->c[k]*xcoefv[j*nc+k]*ycoefv[i*nc+k]; + DFloat y = poly_v->c[0]; for (auto k=1; k< nc; ++k) y+=poly_v->c[k]*coef[j*nc+k]*coef[i*nc+k]; if (y < fl0y) continue; if (y >= flly) continue; SizeT px = x; @@ -392,10 +389,10 @@ template< typename T1, typename T2> for (OMPInt j = 0; j < nRows; ++j) { for (OMPInt i = 0; i < nCols; ++i) { // Compute the original source for this pixel, note order of j and i in P and Q definition of IDL doc. - DFloat x = poly_u->c[0]; for (auto k=1; k< nc; ++k) x+=poly_u->c[k]*xcoefu[j*nc+k]*ycoefu[i*nc+k]; + DFloat x = poly_u->c[0]; for (auto k=1; k< nc; ++k) x+=poly_u->c[k]*coef[j*nc+k]*coef[i*nc+k]; if (x < fl0x) continue; if (x >= fllx) continue; // already initialised to 'missing' value. - DFloat y = poly_v->c[0]; for (auto k=1; k< nc; ++k) y+=poly_v->c[k]*xcoefv[j*nc+k]*ycoefv[i*nc+k]; + DFloat y = poly_v->c[0]; for (auto k=1; k< nc; ++k) y+=poly_v->c[k]*coef[j*nc+k]*coef[i*nc+k]; if (y < fl0y) continue; if (y >= flly) continue; SizeT px = x; @@ -418,8 +415,8 @@ template< typename T1, typename T2> DFloat x = poly_u->c[0]; DFloat y = poly_v->c[0]; for (auto k=1; k< nc; ++k) { - x+=poly_u->c[k]*xcoefu[j*nc+k]*ycoefu[i*nc+k]; - y+=poly_v->c[k]*xcoefv[j*nc+k]*ycoefv[i*nc+k]; + x+=poly_u->c[k]*coef[j*nc+k]*coef[i*nc+k]; + y+=poly_v->c[k]*coef[j*nc+k]*coef[i*nc+k]; } SizeT px=x; SizeT py=y; @@ -437,8 +434,8 @@ template< typename T1, typename T2> DFloat x = poly_u->c[0]; DFloat y = poly_v->c[0]; for (auto k=1; k< nc; ++k) { - x+=poly_u->c[k]*xcoefu[j*nc+k]*ycoefu[i*nc+k]; - y+=poly_v->c[k]*xcoefv[j*nc+k]*ycoefv[i*nc+k]; + x+=poly_u->c[k]*coef[j*nc+k]*coef[i*nc+k]; + y+=poly_v->c[k]*coef[j*nc+k]*coef[i*nc+k]; } SizeT px=x; SizeT py=y; @@ -449,18 +446,7 @@ template< typename T1, typename T2> } } } - free(xcoefu); - free(ycoefu); - free(xcoefv); - free(ycoefv); - free(poly_u->px); - free(poly_u->py); - free(poly_u->c); - free(poly_u); - free(poly_v->px); - free(poly_v->py); - free(poly_v->c); - free(poly_v); + free(coef); return res_; } @@ -469,9 +455,9 @@ template< typename T1, typename T2> const SizeT nCols, const SizeT nRows, BaseGDL* data_, - DFloat * const P, - DFloat * const Q, - DFloat const initvalue_, + const DFloat * P, + const DFloat * Q, + const DFloat initvalue_, const bool doMissing) { // std::cerr << "warp_linear1\n"; SizeT lx = data_->Dim(0); @@ -494,20 +480,20 @@ template< typename T1, typename T2> } } //these are accelerators - not using them increase exec time by 2 or 3. - const DFloat xmax = lx; - const DFloat ymax = ly; - const DFloat xmin = 0; - const DFloat ymin = 0; - const DFloat xbound = lx - 1; //-1 for linear - const DFloat ybound = ly - 1; - const DFloat p0 = P[0]; - const DFloat q0 = Q[0]; - const DFloat p1 = P[1]; - const DFloat q1 = Q[1]; - const DFloat p2 = P[2]; - const DFloat q2 = Q[2]; - const DFloat p3 = P[3]; - const DFloat q3 = Q[3]; + DFloat xmax = lx; + DFloat ymax = ly; + DFloat xmin = 0; + DFloat ymin = 0; + DFloat xbound = lx - 1; //-1 for linear + DFloat ybound = ly - 1; + DFloat p0 = P[0]; + DFloat q0 = Q[0]; + DFloat p1 = P[1]; + DFloat q1 = Q[1]; + DFloat p2 = P[2]; + DFloat q2 = Q[2]; + DFloat p3 = P[3]; + DFloat q3 = Q[3]; DFloat p1j, p3j, q1j, q3j; /* Double loop on the output internal image */ @@ -615,14 +601,14 @@ template< typename T1, typename T2> BaseGDL* warp1( const SizeT nCols, const SizeT nRows, - BaseGDL* const data_, - poly2d* const poly_u, - poly2d* const poly_v, + BaseGDL * data_, + poly2d * poly_u, + poly2d * poly_v, const DFloat initvalue_, - bool doMissing) { + const bool doMissing) { // std::cerr<<"warp1\n"; - const SizeT lx = data_->Dim(0); - const SizeT ly = data_->Dim(1); + SizeT lx = data_->Dim(0); + SizeT ly = data_->Dim(1); dimension dim(nCols, nRows); T1* res_ = new T1(dim, BaseGDL::NOZERO); @@ -646,21 +632,19 @@ template< typename T1, typename T2> DFloat fl0x=0; DFloat fl0y=0; DLong nc=poly_u->nc; - DFloat * const xcoefu=poly2d_compute_init_x(poly_u,lx); - DFloat * const ycoefu=poly2d_compute_init_y(poly_u,lx); - DFloat * const xcoefv=poly2d_compute_init_x(poly_v,ly); - DFloat * const ycoefv=poly2d_compute_init_y(poly_v,ly); + assert(poly_u->nc == poly_v->nc); + DFloat * const coef=poly2d_compute_init(poly_u,lx,ly); //poly_u and poly_va share the same size. /* Double loop on the output image */ if (doMissing) { if ((GDL_NTHREADS = parallelize(nEl)) == 1) { for (OMPInt j = 0; j < nRows; ++j) { for (OMPInt i = 0; i < nCols; ++i) { DFloat x = poly_u->c[0]; - for (auto k = 1; k < nc; ++k) x += poly_u->c[k] * xcoefu[j * nc + k] * ycoefu[i * nc + k]; + for (auto k = 1; k < nc; ++k) x += poly_u->c[k] * coef[j * nc + k] * coef[i * nc + k]; if (x < fl0x) continue; if (x >= fllx) continue; DFloat y = poly_v->c[0]; - for (auto k = 1; k < nc; ++k) y += poly_v->c[k] * xcoefv[j * nc + k] * ycoefv[i * nc + k]; + for (auto k = 1; k < nc; ++k) y += poly_v->c[k] * coef[j * nc + k] * coef[i * nc + k]; if (y < fl0y) continue; if (y >= flly) continue; SizeT px = x; @@ -683,11 +667,11 @@ template< typename T1, typename T2> for (OMPInt j = 0; j < nRows; ++j) { for (OMPInt i = 0; i < nCols; ++i) { DFloat x = poly_u->c[0]; - for (auto k = 1; k < nc; ++k) x += poly_u->c[k] * xcoefu[j * nc + k] * ycoefu[i * nc + k]; + for (auto k = 1; k < nc; ++k) x += poly_u->c[k] * coef[j * nc + k] * coef[i * nc + k]; if (x < fl0x) continue; if (x >= fllx) continue; DFloat y = poly_v->c[0]; - for (auto k = 1; k < nc; ++k) y += poly_v->c[k] * xcoefv[j * nc + k] * ycoefv[i * nc + k]; + for (auto k = 1; k < nc; ++k) y += poly_v->c[k] * coef[j * nc + k] * coef[i * nc + k]; if (y < fl0y) continue; if (y >= flly) continue; SizeT px = x; @@ -718,8 +702,8 @@ template< typename T1, typename T2> DFloat x = poly_u->c[0]; DFloat y = poly_v->c[0]; for (auto k=1; k< nc; ++k) { - x+=poly_u->c[k]*xcoefu[j*nc+k]*ycoefu[i*nc+k]; - y+=poly_v->c[k]*xcoefv[j*nc+k]*ycoefv[i*nc+k]; + x+=poly_u->c[k]*coef[j*nc+k]*coef[i*nc+k]; + y+=poly_v->c[k]*coef[j*nc+k]*coef[i*nc+k]; } SizeT px=x; SizeT py=y; @@ -745,8 +729,8 @@ template< typename T1, typename T2> DFloat x = poly_u->c[0]; DFloat y = poly_v->c[0]; for (auto k=1; k< nc; ++k) { - x+=poly_u->c[k]*xcoefu[j*nc+k]*ycoefu[i*nc+k]; - y+=poly_v->c[k]*xcoefv[j*nc+k]*ycoefv[i*nc+k]; + x+=poly_u->c[k]*coef[j*nc+k]*coef[i*nc+k]; + y+=poly_v->c[k]*coef[j*nc+k]*coef[i*nc+k]; } SizeT px=x; SizeT py=y; @@ -766,18 +750,7 @@ template< typename T1, typename T2> } } } - free(xcoefu); - free(ycoefu); - free(xcoefv); - free(ycoefv); - free(poly_u->px); - free(poly_u->py); - free(poly_u->c); - free(poly_u); - free(poly_v->px); - free(poly_v->py); - free(poly_v->c); - free(poly_v); + free(coef); return res_; } @@ -786,10 +759,10 @@ template< typename T1, typename T2> const SizeT nCols, const SizeT nRows, BaseGDL* data_, - DFloat * const P, - DFloat * const Q, + const DFloat * P, + const DFloat * Q, const DFloat cubicParameter, - DFloat const initvalue_, + const DFloat initvalue_, const bool doMissing) { // std::cerr<<"warp_linear2\n"; const SizeT lx = data_->Dim(0); @@ -841,20 +814,20 @@ template< typename T1, typename T2> leaps[3][3] = 2 + 2 * lx; //these are accelerators - not using them increase exec time by 2 or 3. - const DFloat xmax = lx; - const DFloat ymax = ly; - const DFloat xmin = 1; - const DFloat ymin = 1; - const DFloat xbound = lx - 2; //-2 for cubic - const DFloat ybound = ly - 2; - const DFloat p0 = P[0]; - const DFloat q0 = Q[0]; - const DFloat p1 = P[1]; - const DFloat q1 = Q[1]; - const DFloat p2 = P[2]; - const DFloat q2 = Q[2]; - const DFloat p3 = P[3]; - const DFloat q3 = Q[3]; + DFloat xmax = lx; + DFloat ymax = ly; + DFloat xmin = 1; + DFloat ymin = 1; + DFloat xbound = lx - 2; //-2 for cubic + DFloat ybound = ly - 2; + DFloat p0 = P[0]; + DFloat q0 = Q[0]; + DFloat p1 = P[1]; + DFloat q1 = Q[1]; + DFloat p2 = P[2]; + DFloat q2 = Q[2]; + DFloat p3 = P[3]; + DFloat q3 = Q[3]; DFloat p1j, p3j, q1j, q3j; /* Double loop on the output internal image */ @@ -1001,15 +974,15 @@ template< typename T1, typename T2> BaseGDL* warp2( const SizeT nCols, const SizeT nRows, - BaseGDL* const data_, + BaseGDL * data_, + poly2d * poly_u, + poly2d * poly_v, const DFloat cubicParameter, - poly2d* const poly_u, - poly2d* const poly_v, const DFloat initvalue_, - bool doMissing) { + const bool doMissing) { // std::cerr<<"warp2\n"; - const SizeT lx = data_->Dim(0); - const SizeT ly = data_->Dim(1); + SizeT lx = data_->Dim(0); + SizeT ly = data_->Dim(1); DFloat* kernel=cubicKernel; //already computed bool destroyKernel=false; @@ -1057,31 +1030,28 @@ template< typename T1, typename T2> leaps[3][3] = 2 + 2 * lx; //these are accelerators - not using them increase exec time by 2 or 3. - const DFloat xmax = lx; - const DFloat ymax = ly; - const DFloat xmin = 1; - const DFloat ymin = 1; - const DFloat xbound = lx - 2; //-2 for cubic - const DFloat ybound = ly - 2; + DFloat xmax = lx; + DFloat ymax = ly; + DFloat xmin = 1; + DFloat ymin = 1; + DFloat xbound = lx - 2; //-2 for cubic + DFloat ybound = ly - 2; DLong nc = poly_u->nc; - DFloat * const xcoefu=poly2d_compute_init_x(poly_u,lx); - DFloat * const ycoefu=poly2d_compute_init_y(poly_u,lx); - DFloat * const xcoefv=poly2d_compute_init_x(poly_v,ly); - DFloat * const ycoefv=poly2d_compute_init_y(poly_v,ly); - + assert(poly_u->nc == poly_v->nc); + DFloat * const coef=poly2d_compute_init(poly_u,lx,ly);//poly_u and poly_va share the same size. /* Double loop on the output internal image */ if ((GDL_NTHREADS=parallelize( nEl ))==1) { for (OMPInt j = 0; j < nRows; ++j) { for (OMPInt i = 0; i < nCols; ++i) { // Compute the original source for this pixel, note order of j and i in P and Q definition of IDL doc. DFloat x = poly_u->c[0]; - for (auto k = 1; k < nc; ++k) x += poly_u->c[k] * xcoefu[j * nc + k] * ycoefu[i * nc + k]; + for (auto k = 1; k < nc; ++k) x += poly_u->c[k] * coef[j * nc + k] * coef[i * nc + k]; if (doMissing) { if (x < xmin) continue; if (x >= xmax) continue; } DFloat y = poly_v->c[0]; - for (auto k = 1; k < nc; ++k) y += poly_v->c[k] * xcoefv[j * nc + k] * ycoefv[i * nc + k]; + for (auto k = 1; k < nc; ++k) y += poly_v->c[k] * coef[j * nc + k] * coef[i * nc + k]; if (doMissing) { if (y < ymin) continue; if (y >= ymax) continue; @@ -1141,13 +1111,13 @@ template< typename T1, typename T2> for (OMPInt i = 0; i < nCols; ++i) { // Compute the original source for this pixel, note order of j and i in P and Q definition of IDL doc. DFloat x = poly_u->c[0]; - for (auto k = 1; k < nc; ++k) x += poly_u->c[k] * xcoefu[j * nc + k] * ycoefu[i * nc + k]; + for (auto k = 1; k < nc; ++k) x += poly_u->c[k] * coef[j * nc + k] * coef[i * nc + k]; if (doMissing) { if (x < xmin) continue; if (x >= xmax) continue; } DFloat y = poly_v->c[0]; - for (auto k = 1; k < nc; ++k) y += poly_v->c[k] * xcoefv[j * nc + k] * ycoefv[i * nc + k]; + for (auto k = 1; k < nc; ++k) y += poly_v->c[k] * coef[j * nc + k] * coef[i * nc + k]; if (doMissing) { if (y < ymin) continue; if (y >= ymax) continue; @@ -1201,9 +1171,91 @@ template< typename T1, typename T2> } } } + free(coef); if (destroyKernel) free(kernel); return res_; } +#define POLY2D_FUNCTION warp_linear0 + template< typename T1, typename T2, typename T3, typename T4> // + BaseGDL* poly2d_complex_handle_warp_linear0( + const SizeT nCols, + const SizeT nRows, + T1* data_, + const DFloat * P, + const DFloat * Q, + const DFloat missing, + const bool doMissing) { +#include "snippets/poly2d_complex_handle_function.incpp" + } +#undef POLY2D_FUNCTION +#define POLY2D_FUNCTION warp_linear1 + template< typename T1, typename T2, typename T3, typename T4> // + BaseGDL* poly2d_complex_handle_warp_linear1( + const SizeT nCols, + const SizeT nRows, + T1* data_, + const DFloat * P, + const DFloat * Q, + const DFloat missing, + const bool doMissing) { +#include "snippets/poly2d_complex_handle_function.incpp" + } +#undef POLY2D_FUNCTION +#define POLY2D_FUNCTION warp_linear2 + template< typename T1, typename T2, typename T3, typename T4> // + BaseGDL* poly2d_complex_handle_warp_linear2( + const SizeT nCols, + const SizeT nRows, + T1* data_, + const DFloat * P, + const DFloat * Q, + const DFloat cubicParameter, + const DFloat missing, + const bool doMissing) { +#include "snippets/poly2d_complex_handle_cubic_function.incpp" + } +#undef POLY2D_FUNCTION + // same code for warp, only the P and Q arguments are diffrent type +#define POLY2D_FUNCTION warp0 + template< typename T1, typename T2, typename T3, typename T4> // + BaseGDL* poly2d_complex_handle_warp0( + const SizeT nCols, + const SizeT nRows, + T1* data_, + poly2d * P, + poly2d * Q, + const DFloat missing, + const bool doMissing) { +#include "snippets/poly2d_complex_handle_function.incpp" + } +#undef POLY2D_FUNCTION +#define POLY2D_FUNCTION warp1 + template< typename T1, typename T2, typename T3, typename T4> // + BaseGDL* poly2d_complex_handle_warp1( + const SizeT nCols, + const SizeT nRows, + T1* data_, + poly2d * P, + poly2d * Q, + const DFloat missing, + const bool doMissing) { +#include "snippets/poly2d_complex_handle_function.incpp" + } +#undef POLY2D_FUNCTION +#define POLY2D_FUNCTION warp2 + template< typename T1, typename T2, typename T3, typename T4> // + BaseGDL* poly2d_complex_handle_warp2( + const SizeT nCols, + const SizeT nRows, + T1* data_, + poly2d * P, + poly2d * Q, + const DFloat cubicParameter, + const DFloat missing, + const bool doMissing) { +#include "snippets/poly2d_complex_handle_cubic_function.incpp" + } +#undef POLY2D_FUNCTION BaseGDL* poly_2d_fun(EnvT* e) { @@ -1216,7 +1268,7 @@ template< typename T1, typename T2> BaseGDL* p0 = e->GetParDefined(0); - if (p0->Type() == GDL_COMPLEX || p0->Type() == GDL_COMPLEXDBL) e->Throw("Complex values not supported (FIXME)"); + if (NumericType(p0->Type())==false) e->Throw("Type of operand not supported."); if (p0->Rank() != 2) e->Throw("Array must have 2 dimensions: " + e->GetParString(0)); @@ -1308,6 +1360,12 @@ template< typename T1, typename T2> return warp_linear0< DFloatGDL, DFloat>(nCol, nRow, p0, (DFloat*) P->DataAddr(),(DFloat*) Q->DataAddr(), missing, doMissing); } else if (p0->Type() == GDL_DOUBLE) { return warp_linear0< DDoubleGDL, DDouble>(nCol, nRow, p0, (DFloat*) P->DataAddr(),(DFloat*) Q->DataAddr(), missing, doMissing); + } else if (p0->Type() == GDL_COMPLEX) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexGDL* p0c=static_cast(p0); + return poly2d_complex_handle_warp_linear0(nCol, nRow, p0c, (DFloat*) P->DataAddr(),(DFloat*) Q->DataAddr(), missing, doMissing); + } else if (p0->Type() == GDL_COMPLEXDBL) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexDblGDL* p0c=static_cast(p0); + return poly2d_complex_handle_warp_linear0(nCol, nRow, p0c, (DFloat*) P->DataAddr(),(DFloat*) Q->DataAddr(), missing, doMissing); } } else if (interp==1) { if (p0->Type() == GDL_BYTE) { @@ -1328,7 +1386,13 @@ template< typename T1, typename T2> return warp_linear1< DFloatGDL, DFloat>(nCol, nRow, p0, (DFloat*) P->DataAddr(),(DFloat*) Q->DataAddr(), missing, doMissing); } else if (p0->Type() == GDL_DOUBLE) { return warp_linear1< DDoubleGDL, DDouble>(nCol, nRow, p0, (DFloat*) P->DataAddr(),(DFloat*) Q->DataAddr(), missing, doMissing); - } + } else if (p0->Type() == GDL_COMPLEX) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexGDL* p0c = static_cast (p0); + return poly2d_complex_handle_warp_linear1(nCol, nRow, p0c, (DFloat*) P->DataAddr(), (DFloat*) Q->DataAddr(), missing, doMissing); + } else if (p0->Type() == GDL_COMPLEXDBL) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexDblGDL* p0c = static_cast (p0); + return poly2d_complex_handle_warp_linear1(nCol, nRow, p0c, (DFloat*) P->DataAddr(), (DFloat*) Q->DataAddr(), missing, doMissing); + } } else if (interp==2) { if (p0->Type() == GDL_BYTE) { return warp_linear2< DByteGDL, DByte>(nCol, nRow, p0, (DFloat*) P->DataAddr(),(DFloat*) Q->DataAddr(), cubicParameter, missing, doMissing); @@ -1348,103 +1412,136 @@ template< typename T1, typename T2> return warp_linear2< DFloatGDL, DFloat>(nCol, nRow, p0, (DFloat*) P->DataAddr(),(DFloat*) Q->DataAddr(), cubicParameter, missing, doMissing); } else if (p0->Type() == GDL_DOUBLE) { return warp_linear2< DDoubleGDL, DDouble>(nCol, nRow, p0, (DFloat*) P->DataAddr(),(DFloat*) Q->DataAddr(), cubicParameter, missing, doMissing); - } + } else if (p0->Type() == GDL_COMPLEX) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexGDL* p0c = static_cast (p0); + return poly2d_complex_handle_warp_linear2(nCol, nRow, p0c, (DFloat*) P->DataAddr(), (DFloat*) Q->DataAddr(), cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_COMPLEXDBL) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexDblGDL* p0c = static_cast (p0); + return poly2d_complex_handle_warp_linear2(nCol, nRow, p0c, (DFloat*) P->DataAddr(), (DFloat*) Q->DataAddr(), cubicParameter, missing, doMissing); + } } } - + else { - - //NON-LINEAR Polynomial - poly2d* poly_u; - poly2d* poly_v; + //NON-LINEAR Polynomial - poly_u = (poly2d *) malloc(sizeof (poly2d)); - poly_u->nc = nc; - poly_u->px = (DLong *) malloc(nc * sizeof (DLong)); - poly_u->py = (DLong *) malloc(nc * sizeof (DLong)); - poly_u->c = (DFloat *) malloc(nc * sizeof (DFloat)); + poly2d* poly_u; + poly2d* poly_v; - for (SizeT i = 0; i < (nDegree + 1)*(nDegree + 1); ++i) { - poly_u->px[i] = i / (nDegree + 1); - poly_u->py[i] = i - (poly_u->px[i] * (nDegree + 1)); - poly_u->c[i] = (*P)[poly_u->px[i]+(nDegree + 1) * poly_u->py[i]]; - } + poly_u = (poly2d *) malloc(sizeof (poly2d)); + poly_u->nc = nc; + poly_u->px = (DLong *) malloc(nc * sizeof (DLong)); + poly_u->py = (DLong *) malloc(nc * sizeof (DLong)); + poly_u->c = (DFloat *) malloc(nc * sizeof (DFloat)); - poly_v = (poly2d *) malloc(sizeof (poly2d)); - poly_v->nc = nc; - poly_v->px = (DLong *) malloc(nc * sizeof (DLong)); - poly_v->py = (DLong *) malloc(nc * sizeof (DLong)); - poly_v->c = (DFloat *) malloc(nc * sizeof (DFloat)); + for (SizeT i = 0; i < (nDegree + 1)*(nDegree + 1); ++i) { + poly_u->px[i] = i / (nDegree + 1); + poly_u->py[i] = i - (poly_u->px[i] * (nDegree + 1)); + poly_u->c[i] = (*P)[poly_u->px[i]+(nDegree + 1) * poly_u->py[i]]; + } - for (SizeT i = 0; i < (nDegree + 1)*(nDegree + 1); ++i) { - poly_v->px[i] = i / (nDegree + 1); - poly_v->py[i] = i - (poly_v->px[i] * (nDegree + 1)); - poly_v->c[i] = (*Q)[poly_v->px[i]+(nDegree + 1) * poly_v->py[i]]; - } - if (interp==0) { - if (p0->Type() == GDL_BYTE) { - return warp0< DByteGDL, DByte>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_INT) { - return warp0< DIntGDL, DInt>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_UINT) { - return warp0< DUIntGDL, DUInt>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_LONG) { - return warp0< DLongGDL, DLong>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_ULONG) { - return warp0< DULongGDL, DULong>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_LONG64) { - return warp0< DLong64GDL, DLong64>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_ULONG64) { - return warp0< DULong64GDL, DULong64>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_FLOAT) { - return warp0< DFloatGDL, DFloat>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_DOUBLE) { - return warp0< DDoubleGDL, DDouble>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } - } else if (interp==1) { - if (p0->Type() == GDL_BYTE) { - return warp1< DByteGDL, DByte>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_INT) { - return warp1< DIntGDL, DInt>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_UINT) { - return warp1< DUIntGDL, DUInt>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_LONG) { - return warp1< DLongGDL, DLong>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_ULONG) { - return warp1< DULongGDL, DULong>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_LONG64) { - return warp1< DLong64GDL, DLong64>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_ULONG64) { - return warp1< DULong64GDL, DULong64>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_FLOAT) { - return warp1< DFloatGDL, DFloat>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_DOUBLE) { - return warp1< DDoubleGDL, DDouble>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); - } - } else if (interp==2) { - if (p0->Type() == GDL_BYTE) { - return warp2< DByteGDL, DByte>(nCol, nRow, p0, cubicParameter, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_INT) { - return warp2< DIntGDL, DInt>(nCol, nRow, p0, cubicParameter, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_UINT) { - return warp2< DUIntGDL, DUInt>(nCol, nRow, p0, cubicParameter, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_LONG) { - return warp2< DLongGDL, DLong>(nCol, nRow, p0, cubicParameter, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_ULONG) { - return warp2< DULongGDL, DULong>(nCol, nRow, p0, cubicParameter, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_LONG64) { - return warp2< DLong64GDL, DLong64>(nCol, nRow, p0, cubicParameter, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_ULONG64) { - return warp2< DULong64GDL, DULong64>(nCol, nRow, p0, cubicParameter, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_FLOAT) { - return warp2< DFloatGDL, DFloat>(nCol, nRow, p0, cubicParameter, poly_u, poly_v, missing, doMissing); - } else if (p0->Type() == GDL_DOUBLE) { - return warp2< DDoubleGDL, DDouble>(nCol, nRow, p0, cubicParameter, poly_u, poly_v, missing, doMissing); - } - } + poly_v = (poly2d *) malloc(sizeof (poly2d)); + poly_v->nc = nc; + poly_v->px = (DLong *) malloc(nc * sizeof (DLong)); + poly_v->py = (DLong *) malloc(nc * sizeof (DLong)); + poly_v->c = (DFloat *) malloc(nc * sizeof (DFloat)); - } + for (SizeT i = 0; i < (nDegree + 1)*(nDegree + 1); ++i) { + poly_v->px[i] = i / (nDegree + 1); + poly_v->py[i] = i - (poly_v->px[i] * (nDegree + 1)); + poly_v->c[i] = (*Q)[poly_v->px[i]+(nDegree + 1) * poly_v->py[i]]; + } + BaseGDL* ret; + if (interp == 0) { + if (p0->Type() == GDL_BYTE) { + ret = warp0< DByteGDL, DByte>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_INT) { + ret = warp0< DIntGDL, DInt>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_UINT) { + ret = warp0< DUIntGDL, DUInt>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_LONG) { + ret = warp0< DLongGDL, DLong>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_ULONG) { + ret = warp0< DULongGDL, DULong>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_LONG64) { + ret = warp0< DLong64GDL, DLong64>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_ULONG64) { + ret = warp0< DULong64GDL, DULong64>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_FLOAT) { + ret = warp0< DFloatGDL, DFloat>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_DOUBLE) { + ret = warp0< DDoubleGDL, DDouble>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_COMPLEX) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexGDL* p0c = static_cast (p0); + ret = poly2d_complex_handle_warp0(nCol, nRow, p0c, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_COMPLEXDBL) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexDblGDL* p0c = static_cast (p0); + ret = poly2d_complex_handle_warp0(nCol, nRow, p0c, poly_u, poly_v, missing, doMissing); + } + } else if (interp == 1) { + if (p0->Type() == GDL_BYTE) { + ret = warp1< DByteGDL, DByte>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_INT) { + ret = warp1< DIntGDL, DInt>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_UINT) { + ret = warp1< DUIntGDL, DUInt>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_LONG) { + ret = warp1< DLongGDL, DLong>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_ULONG) { + ret = warp1< DULongGDL, DULong>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_LONG64) { + ret = warp1< DLong64GDL, DLong64>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_ULONG64) { + ret = warp1< DULong64GDL, DULong64>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_FLOAT) { + ret = warp1< DFloatGDL, DFloat>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_DOUBLE) { + ret = warp1< DDoubleGDL, DDouble>(nCol, nRow, p0, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_COMPLEX) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexGDL* p0c = static_cast (p0); + ret = poly2d_complex_handle_warp1(nCol, nRow, p0c, poly_u, poly_v, missing, doMissing); + } else if (p0->Type() == GDL_COMPLEXDBL) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexDblGDL* p0c = static_cast (p0); + ret = poly2d_complex_handle_warp1(nCol, nRow, p0c, poly_u, poly_v, missing, doMissing); + } + } else if (interp == 2) { + if (p0->Type() == GDL_BYTE) { + ret = warp2< DByteGDL, DByte>(nCol, nRow, p0, poly_u, poly_v, cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_INT) { + ret = warp2< DIntGDL, DInt>(nCol, nRow, p0, poly_u, poly_v, cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_UINT) { + ret = warp2< DUIntGDL, DUInt>(nCol, nRow, p0, poly_u, poly_v, cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_LONG) { + ret = warp2< DLongGDL, DLong>(nCol, nRow, p0, poly_u, poly_v, cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_ULONG) { + ret = warp2< DULongGDL, DULong>(nCol, nRow, p0, poly_u, poly_v, cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_LONG64) { + ret = warp2< DLong64GDL, DLong64>(nCol, nRow, p0, poly_u, poly_v, cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_ULONG64) { + ret = warp2< DULong64GDL, DULong64>(nCol, nRow, p0, poly_u, poly_v, cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_FLOAT) { + ret = warp2< DFloatGDL, DFloat>(nCol, nRow, p0, poly_u, poly_v, cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_DOUBLE) { + ret = warp2< DDoubleGDL, DDouble>(nCol, nRow, p0, poly_u, poly_v, cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_COMPLEX) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexGDL* p0c = static_cast (p0); + ret = poly2d_complex_handle_warp2(nCol, nRow, p0c, poly_u, poly_v, cubicParameter, missing, doMissing); + } else if (p0->Type() == GDL_COMPLEXDBL) { //duplicating warp code for Complex would make the whole thing unmaintainable. Use a simpler solution. + DComplexDblGDL* p0c = static_cast (p0); + ret = poly2d_complex_handle_warp2(nCol, nRow, p0c, poly_u, poly_v, cubicParameter, missing, doMissing); + } + } + free(poly_u->px); + free(poly_u->py); + free(poly_u->c); + free(poly_u); + free(poly_v->px); + free(poly_v->py); + free(poly_v->c); + free(poly_v); + return ret; + } return NULL; } } diff --git a/src/snippets/poly2d_complex_handle_cubic_function.incpp b/src/snippets/poly2d_complex_handle_cubic_function.incpp new file mode 100644 index 000000000..d5fa79011 --- /dev/null +++ b/src/snippets/poly2d_complex_handle_cubic_function.incpp @@ -0,0 +1,18 @@ + T2* tmp = new T2(data_->Dim(), BaseGDL::NOZERO); + for (auto i = 0; i < data_->N_Elements(); ++i) (*tmp)[i] = (*data_)[i].real(); + BaseGDL* rea = POLY2D_FUNCTION(nCols, nRows, tmp, P, Q, cubicParameter, missing, doMissing); + for (auto i = 0; i < data_->N_Elements(); ++i) (*tmp)[i] = (*data_)[i].imag(); + BaseGDL* ima = POLY2D_FUNCTION(nCols, nRows, tmp, P, Q, cubicParameter, 0, doMissing); //missing imaginary is always 0 in IDL. + GDLDelete(tmp); + dimension dim(nCols, nRows); + T1* comp = new T1(dim, BaseGDL::NOZERO); + T3* fcomp = static_cast (comp->DataAddr()); + T4* realval = static_cast (rea->DataAddr()); + T4* imagval = static_cast (ima->DataAddr()); + for (auto i = 0; i < nCols * nRows; ++i) { + reinterpret_cast(fcomp)[2 * i] = realval[i]; + reinterpret_cast(fcomp)[2 * i +1] = imagval[i]; + } + GDLDelete(rea); + GDLDelete(ima); + return comp; diff --git a/src/snippets/poly2d_complex_handle_function.incpp b/src/snippets/poly2d_complex_handle_function.incpp new file mode 100644 index 000000000..c401128ca --- /dev/null +++ b/src/snippets/poly2d_complex_handle_function.incpp @@ -0,0 +1,18 @@ + T2* tmp = new T2(data_->Dim(), BaseGDL::NOZERO); + for (auto i = 0; i < data_->N_Elements(); ++i) (*tmp)[i] = (*data_)[i].real(); + BaseGDL* rea = POLY2D_FUNCTION(nCols, nRows, tmp, P, Q, missing, doMissing); + for (auto i = 0; i < data_->N_Elements(); ++i) (*tmp)[i] = (*data_)[i].imag(); + BaseGDL* ima = POLY2D_FUNCTION(nCols, nRows, tmp, P, Q, 0, doMissing); //missing imaginary is always 0 in IDL. + GDLDelete(tmp); + dimension dim(nCols, nRows); + T1* comp = new T1(dim, BaseGDL::NOZERO); + T3* fcomp = static_cast (comp->DataAddr()); + T4* realval = static_cast (rea->DataAddr()); + T4* imagval = static_cast (ima->DataAddr()); + for (auto i = 0; i < nCols * nRows; ++i) { + reinterpret_cast(fcomp)[2 * i] = realval[i]; + reinterpret_cast(fcomp)[2 * i +1] = imagval[i]; + } + GDLDelete(rea); + GDLDelete(ima); + return comp; diff --git a/testsuite/test_all_basic_functions.pro b/testsuite/test_all_basic_functions.pro index b0abc3fad..f20b1b99f 100644 --- a/testsuite/test_all_basic_functions.pro +++ b/testsuite/test_all_basic_functions.pro @@ -163,7 +163,7 @@ if (section eq 0 or section eq 5) then begin 'print,what[i] & for k=0,all_numeric do begin & subclock=tic(typenames[k]) & ret=atan(*big[k],*big[k]) & toc,subclock & end ',$ 'print,what[i] & kernel=[ [0,1,0],[-1,0,1],[0,-1,0] ] & for k=0,all_numeric do begin & subclock=tic(typenames[k]) & ret=convol(*big[k],kernel) & toc,subclock & end ',$ 'print,what[i] & z=findgen(size) & for k=0,all_numeric do begin & subclock=tic(typenames[k]) & ret=interpolate(*big[k],z) & toc,subclock & end ',$ - 'print,what[i] & for k=0,not_complex do begin & subclock=tic(typenames[k]) & ret=poly_2d(*big[k],p,q) & toc,subclock & end ',$ + 'print,what[i] & for k=0,all_numeric do begin & subclock=tic(typenames[k]) & ret=poly_2d(*big[k],p,q) & toc,subclock & end ',$ 'print,what[i] & for k=0,all_numeric do begin & subclock=tic(typenames[k]) & tvscl,(*big[k]) & toc,subclock & end ' ] for i=0,n_elements(calls)-1 do begin & clock=tic(what[i]) & z=execute(calls[i]) & toc,clock & end diff --git a/testsuite/test_bugs_poly2d.pro b/testsuite/test_bugs_poly2d.pro index aca7d7833..5482f55e3 100644 --- a/testsuite/test_bugs_poly2d.pro +++ b/testsuite/test_bugs_poly2d.pro @@ -5,7 +5,7 @@ if KEYWORD_SET(help) then begin print, ' no_exit=no_exit, test=test' return endif -; +; checking errors for interpolation 0 ONLY - fixme. nb_errors=0 a=bindgen(3,3) degree=0 @@ -59,6 +59,41 @@ if KEYWORD_SET(test) then STOP,nb_errors ; BANNER_FOR_TESTSUITE, 'TEST_BUGS_POLY2D', nb_errors ; +; ADDING unchecked commands to serve as a coverage test for all cases: +; initialisations: floats at end, since some commands do not accpet floats/doubles/complex + typecodes=[1,2,3,12,13,14,15,4,5,6,9] + typenames=[" BYTE"," INT"," LONG"," UINT"," ULONG"," LONG64"," ULONG64"," FLOAT"," DOUBLE"," COMPLEX"," DCOMPLEX"] + all_numeric=10 + seed=33 + ; use non-trivial dimensions + dim1=7 + dim2=9 + a=randomn(seed,dim1,dim2,/double)*randomu(seed,dim1,dim2,/ulong) + big=ptrarr(11,/allo) + k=0 & foreach i,typecodes do begin & *big[k]=fix(a,type=i) & k++ &end +; make dummy poly_2D on all types: 1) dimension 1 +; the linear case without the x*y coefficient is already tested by test_all_basic_functions. +p=replicate(0.3,4) & q=replicate(0.4,4) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,0, 27, 19) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,1, 27, 19) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,2, 27, 19) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,2, 27, 19, cub=-0.5) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,0, 27, 19, miss=-1) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,1, 27, 19, miss=-1) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,2, 27, 19, miss=-1) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,2, 27, 19, cub=-0.5, miss=-1) +; idem dimension 2 +p=replicate(0.3,9) & q=replicate(0.4,9) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,0, 27, 19) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,1, 27, 19) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,2, 27, 19) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,2, 27, 19, cub=-0.5) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,0, 27, 19, miss=-1) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,1, 27, 19, miss=-1) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,2, 27, 19, miss=-1) + for k=0,all_numeric do ret=poly_2d(*big[k],p,q,2, 27, 19, cub=-0.5, miss=-1) + + if (nb_errors GT 0) AND ~KEYWORD_SET(no_exit) then EXIT, status=1 ; if KEYWORD_SET(test) then STOP