From 1d03fa5f4a8b4ce33ea655c2146b0f46dfb946c6 Mon Sep 17 00:00:00 2001 From: oguyon Date: Sat, 18 Nov 2023 23:41:39 -1000 Subject: [PATCH] added image unfolding --- .../linalgebra/basis_rotate_match.c | 261 ++++++++------- src/COREMOD_arith/CMakeLists.txt | 2 + src/COREMOD_arith/COREMOD_arith.c | 3 + src/COREMOD_arith/image_unfold.c | 312 ++++++++++++++++++ src/COREMOD_arith/image_unfold.h | 6 + 5 files changed, 452 insertions(+), 132 deletions(-) create mode 100644 src/COREMOD_arith/image_unfold.c create mode 100644 src/COREMOD_arith/image_unfold.h diff --git a/plugins/milk-extra-src/linalgebra/basis_rotate_match.c b/plugins/milk-extra-src/linalgebra/basis_rotate_match.c index f943b908..308fe203 100644 --- a/plugins/milk-extra-src/linalgebra/basis_rotate_match.c +++ b/plugins/milk-extra-src/linalgebra/basis_rotate_match.c @@ -98,12 +98,12 @@ errno_t compute_basis_rotate_match( // internal Arot array, double for improved precision // - double * __restrict Arot = (double *) malloc(sizeof(double) * Adim * Adim ); + double * Arot = (double *) malloc(sizeof(double) * Adim * Adim ); // internal copy of imginAB, double for improved precision // - double * __restrict matAB = (double *) malloc(sizeof(double) * Adim * Bdim); + double * matAB = (double *) malloc(sizeof(double) * Adim * Bdim); // loop stop condition: toggles to 0 when done @@ -111,7 +111,7 @@ errno_t compute_basis_rotate_match( // diagonal values // - double * __restrict diagVal = (double*) malloc(sizeof(double)*Adim); + double * diagVal = (double*) malloc(sizeof(double)*Adim); double diagVal_lim = 0.0; int loopiter = 0; @@ -133,6 +133,7 @@ errno_t compute_basis_rotate_match( while( loopOK == 1 ) { + // Initialize: set Arot to identity matrix // for(uint64_t ii=0; ii 0.0 ) + { + dcoeff *= posSideAmp; + } + optall += dcoeff * matAB[iia*Bdim + iib] * dcoeff * matAB[iia*Bdim + iib]; } } - printf("iter %4d val = %g\n", loopiter, optall); + printf("iter %4d / %4d dangle = %f val = %g\n", loopiter, loopiterMax, dangle, optall); - for ( int n0 = 0; n0 < Adim; n0++) //Adim; n0++ ) + for ( int n0 = 0; n0 < Adim; n0++) { for ( int n1 = n0+1; n1 < Adim; n1++ ) { // testing rotation n0 n1, dangle // ref value - double optvalN0 = 0.0; // to be minimized - double optvalD0 = 0.0; // to be maximized - double optval0 = 0.0; // equal o optvalN/optvalD - ratio to be minimized - - double optvalN1 = 0.0; // to be minimized - double optvalD1 = 0.0; // to be maximized - double optval1 = 0.0; // equal o optvalN/optvalD - ratio to be minimized + double optval0 = 0.0; // to be minimized + double optval1 = 0.0; // to be minimized + double optvalpos0 = 0.0; + double optvalpos1 = 0.0; - - double optvalpos0 = 0.0; - double optvalNpos0 = 0.0; - double optvalDpos0 = 0.0; - - double optvalpos1 = 0.0; - double optvalNpos1 = 0.0; - double optvalDpos1 = 0.0; - - - - double optvalneg0 = 0.0; - double optvalNneg0 = 0.0; - double optvalDneg0 = 0.0; - - double optvalneg1 = 0.0; - double optvalNneg1 = 0.0; - double optvalDneg1 = 0.0; + double optvalneg0 = 0.0; + double optvalneg1 = 0.0; for(uint32_t ii=0; ii n0 ) - if ( dx0 > 0 ) - { - dcoeff0 = pow(dx0, alphap); - //dcoeff0 = ii-n0; - } - if(dcoeff0 > alphalim) - { - dcoeff0 = alphalim; - } + double dcoeff0 = pow(dx0*dx0, alphap); + double dcoeff1 = pow(dx1*dx1, alphap); - double dcoeff1 = 0.0; - //if ( ii > n1 ) - if ( dx1 > 0 ) + if( dx0 > 0.0 ) { - dcoeff1 = pow(dx1, alphap); - //dcoeff1 = ii-n1; + dcoeff0 *= posSideAmp; } - if(dcoeff1 > alphalim) + if( dx1 > 0.0 ) { - dcoeff1 = alphalim; + dcoeff1 *= posSideAmp; } - optvalN0 += dcoeff0 * v0*v0; - optvalN1 += dcoeff1 * v1*v1; - - // maximize sum squared of coefficitents to the left (and including) diagonal - // this ensures the target modes are represented - //if(ii <= n0) - if ( dx0 < 0.0 ) - { - optvalD0 += v0*v0; - //printf(" %4d optval1D += %g (%g) -> %g\n", ii, v0*v0, v0, optvalD); - } - //if(ii <= n1) - if ( dx1 < 0.0 ) - { - optvalD1 += v1*v1; - //printf(" %4d optval1D += %g (%g) -> %g\n", ii, v1*v1, v1, optvalD); - } + double v0 = matAB[n0*Bdim + ii]; + double v1 = matAB[n1*Bdim + ii]; + // optimization metric without rotation + optval0 += dcoeff0 * v0*v0; + optval1 += dcoeff1 * v1*v1; // perform rotation, weite to n0array and n1array n0arraypos[ii] = v0 * cos(dangle) - v1 * sin(dangle); @@ -448,47 +415,18 @@ errno_t compute_basis_rotate_match( n0arrayneg[ii] = v0 * cos(dangle) + v1 * sin(dangle); n1arrayneg[ii] = - v0 * sin(dangle) + v1 * cos(dangle); - optvalNpos0 += dcoeff0 * n0arraypos[ii] * n0arraypos[ii]; - optvalNpos1 += dcoeff1 * n1arraypos[ii] * n1arraypos[ii]; - if ( dx0 < 0.0 ) - { - optvalDpos0 += n0arraypos[ii] * n0arraypos[ii]; - } - if ( dx1 < 0.0 ) - { - optvalDpos1 += n1arraypos[ii] * n1arraypos[ii]; - } + // optimization metric with positive rotation + optvalpos0 += dcoeff0 * n0arraypos[ii] * n0arraypos[ii]; + optvalpos1 += dcoeff1 * n1arraypos[ii] * n1arraypos[ii]; + // optimization metric with negative rotation + optvalneg0 += dcoeff0 * n0arrayneg[ii] * n0arrayneg[ii]; + optvalneg1 += dcoeff1 * n1arrayneg[ii] * n1arrayneg[ii]; - optvalNneg0 += dcoeff0 * n0arrayneg[ii] * n0arrayneg[ii]; - optvalNneg1 += dcoeff1 * n1arrayneg[ii] * n1arrayneg[ii]; - if ( dx0 < 0.0 ) - { - optvalDneg0 += n0arrayneg[ii] * n0arrayneg[ii]; - } - if ( dx1 < 0.0 ) - { - optvalDneg1 += n1arrayneg[ii] * n1arrayneg[ii]; - } } - double epsN = 1e-8; - double epsD = 1e-16;// avoid division by zero - - optval0 = (optvalN0 + epsN) / (optvalD0 + epsD); - optvalpos0 = (optvalNpos0 + epsN) / (optvalDpos0 + epsD); - optvalneg0 = (optvalNneg0 + epsN) / (optvalDneg0 + epsD); - - optval1 = (optvalN1 + epsN) / (optvalD1 + epsD); - optvalpos1 = (optvalNpos1 + epsN) / (optvalDpos1 + epsD); - optvalneg1 = (optvalNneg1 + epsN) / (optvalDneg1 + epsD); - - // quantity to be minimized - double n0coeff = 1.0; ///pow(1.0+n0, 2.0); - double n1coeff = 1.0; ///pow(1.0+n1, 2.0); - - double optval = optval0*n0coeff;// + optval1*n1coeff; - double optvalneg = optvalneg0*n0coeff;// + optvalneg1*n1coeff; - double optvalpos = optvalpos0*n0coeff;// + optvalpos1*n1coeff; + double optval = optval0 + optval1; + double optvalneg = optvalneg0 + optvalneg1; + double optvalpos = optvalpos0 + optvalpos1; //printf(" [%3d - %3d] %g %g %g\n", n0, n1, optvalneg, optval, optvalpos); @@ -502,7 +440,7 @@ errno_t compute_basis_rotate_match( } else if(optvalpos < optval) { - // rotate neg + // rotate pos optrotangle = dangle; cntpos++; } @@ -565,15 +503,74 @@ errno_t compute_basis_rotate_match( } printf(" [%5ld %5ld %5ld] %g\n", cntneg, cntmid, cntpos, dangle); - if( cntmid > 30*(cntpos+cntneg) ) + if( cntmid > 100.0*(cntpos+cntneg) ) + { + dangle *= 0.98; + } + + + + + + // Measure, for each A mode, the effective index of B modes + // + long * iarray = (long *) malloc(sizeof(long) * Adim); + for( int iia = 0; iia < Adim; iia++ ) + { + iarray[iia] = iia; + double Beff = 0.0; + double Beffcnt = 0.0; + for(uint32_t iib=0; iibnaxis = 2; imgArot->size[0] = Adim; imgArot->size[1] = Adim; @@ -599,7 +596,7 @@ errno_t compute_basis_rotate_match( free(Arot); - //copy matAB to ouput +//copy matAB to ouput for(uint64_t ii = 0; iiarray.F[ii] = matAB[ii]; diff --git a/src/COREMOD_arith/CMakeLists.txt b/src/COREMOD_arith/CMakeLists.txt index 35b029eb..0ed72d09 100644 --- a/src/COREMOD_arith/CMakeLists.txt +++ b/src/COREMOD_arith/CMakeLists.txt @@ -21,6 +21,7 @@ set(SOURCEFILES image_slicenormalize.c image_stats.c image_total.c + image_unfold.c image_dxdy.c imfunctions.c mathfuncs.c @@ -48,6 +49,7 @@ set(INCLUDEFILES image_slicenormalize.h image_stats.h image_total.h + image_unfold.h image_dxdy.h imfunctions.h mathfuncs.h diff --git a/src/COREMOD_arith/COREMOD_arith.c b/src/COREMOD_arith/COREMOD_arith.c index b8ec5107..40b6a84d 100644 --- a/src/COREMOD_arith/COREMOD_arith.c +++ b/src/COREMOD_arith/COREMOD_arith.c @@ -62,6 +62,7 @@ #include "image_set_col.h" #include "image_set_row.h" #include "image_setzero.h" +#include "image_unfold.h" #include "image_pixremap.h" #include "image_pixunmap.h" @@ -111,6 +112,8 @@ static errno_t init_module_CLI() CLIADDCMD_COREMODE_arith__pixremap(); CLIADDCMD_COREMODE_arith__pixunmap(); + CLIADDCMD_COREMOD_arith__image_unfold(); + return RETURN_SUCCESS; } diff --git a/src/COREMOD_arith/image_unfold.c b/src/COREMOD_arith/image_unfold.c new file mode 100644 index 00000000..b89b76b3 --- /dev/null +++ b/src/COREMOD_arith/image_unfold.c @@ -0,0 +1,312 @@ +#include + +#include "CommandLineInterface/CLIcore.h" +#include "COREMOD_memory/COREMOD_memory.h" + + + +// input image names +static char *inimname; +static char *maskimname; + +static char *outimname; + + +static uint32_t *axisA; +static long fpi_axisA = -1; + +static uint32_t *axisB; +static long fpi_axisB = -1; + +static uint32_t *colsize; +static long fpi_colsize = -1; + +static char *auxin; + +static uint64_t *modeRMS; + + + +static CLICMDARGDEF farg[] = +{ + { + CLIARG_IMG, + ".inim", + "input image", + "im0", + CLIARG_VISIBLE_DEFAULT, + (void **) &inimname, + NULL + }, + { + CLIARG_STR, + ".outim", + "output image", + "imout", + CLIARG_VISIBLE_DEFAULT, + (void **) &outimname, + NULL + }, + { + CLIARG_UINT32, + ".axisA", + "axis to merged", + "2", + CLIARG_VISIBLE_DEFAULT, + (void **) &axisA, + &fpi_axisA + }, + { + CLIARG_UINT32, + ".axisB", + "merge into this axis", + "0", + CLIARG_VISIBLE_DEFAULT, + (void **) &axisB, + &fpi_axisB + }, + { + CLIARG_UINT32, + ".colsize", + "column size", + "10", + CLIARG_HIDDEN_DEFAULT, + (void **) &colsize, + &fpi_colsize + } +}; + + + + +static CLICMDDATA CLIcmddata = +{ + "unfold", + "image unfold, merge axis A into axis B", + CLICMD_FIELDS_DEFAULTS +}; + + + +// detailed help +static errno_t help_function() +{ + printf("Unfold image: redude number of axis\n"); + printf("Example, 3D image input [X,Y,Z]\n"); + printf("merge axis 2 into axis 0 -> size [XZ, Y]\n"); + printf("if colsize specified, array size is [X*colsize, Y*n]\n"); + printf("with n sufficiently large to include all pixels\n"); + + return RETURN_SUCCESS; +} + + + + +errno_t image_unfold( + IMGID inimg, + IMGID *outimg, + uint8_t axisA, + uint8_t axisB, + int colsize +) +{ + DEBUG_TRACE_FSTART(); + + resolveIMGID(&inimg, ERRMODE_ABORT); + + + resolveIMGID(outimg, ERRMODE_NULL); + if( outimg->ID == -1) + { + copyIMGID(&inimg, outimg); + } + + // output image size + outimg->naxis = inimg.md->naxis - 1; + + + // remove missing axis + { + uint8_t axout = 0; + for( uint8_t axin=0; axinnaxis; axin++) + { + if( axin != axisA ) + { + outimg->size[axout] = inimg.md->size[axin]; + axout ++; + } + } + } + + // destination axis to grow + uint8_t axis0 = 0; + if( axisA > axisB ) + { + axis0 = axisB; + } + else + { + axis0 = axisB-1; + } + + // overflow destination axis to grow + uint8_t axis1 = 0; + uint8_t axisC = 0; // in input image + if( (axis0 == 0 ) && (outimg->naxis >1) ) + { + axis1 = 1; + axisC = 1; + } + + + int mdimsize = 0; + if( axis0 == axis1 ) + { + outimg->size[axis0] *= inimg.md->size[axisA]; + } + else + { + int mdim0 = 0; // multiplicative on axis0 + int mdim1 = 1; // multiplicative on axis1 + + for( uint32_t ii=0; iisize[axisA]; ii++) + { + mdim0 ++; + if(mdim0 == colsize) + { + mdim0 = 0; + mdim1++; + } + } + if(mdim1 > 0) + { + mdim0 = colsize; + } + + outimg->size[axis0] *= mdim0; + outimg->size[axis1] *= mdim1; + + mdimsize = inimg.md->size[axisC] * outimg->size[axis0]; + } + + createimagefromIMGID(outimg); + + + + + // copy data to ouput + + // destination pix coord + uint32_t ii = 0; + uint32_t jj = 0; + + uint64_t pixi = 0; + uint64_t pixo = 0; + for( uint32_t pixi2=0; pixi2 < inimg.md->size[2]; pixi2++) + { + for( uint32_t pixi1=0; pixi1 < inimg.md->size[1]; pixi1++) + { + for( uint32_t pixi0=0; pixi0 < inimg.md->size[0]; pixi0++) + { + pixo = jj; + pixo *= outimg->md->size[0]; + + pixo += ii % outimg->md->size[0]; + pixo += mdimsize * ( ii / outimg->md->size[0] ); + + + outimg->im->array.F[pixo] = inimg.im->array.F[pixi]; + pixi ++; + + ii++; + } + if (( axisA == 1) && ( axisB == 0)) + { + // do nothing + } + else + { + jj ++; + ii -= inimg.md->size[0]; + } + } + + if (( axisA == 2) && ( axisB == 1)) + { + // do nothing + } + else + { + ii += inimg.md->size[0]; + jj -= inimg.md->size[1]; + } + + + } + + + + + + + + + + DEBUG_TRACE_FEXIT(); + return RETURN_SUCCESS; +} + + + + + + +static errno_t compute_function() +{ + DEBUG_TRACE_FSTART(); + + IMGID inimg = mkIMGID_from_name(inimname); + resolveIMGID(&inimg, ERRMODE_ABORT); + + IMGID outimg = mkIMGID_from_name(outimname); + + INSERT_STD_PROCINFO_COMPUTEFUNC_INIT + + + INSERT_STD_PROCINFO_COMPUTEFUNC_LOOPSTART + { + + image_unfold( + inimg, + &outimg, + *axisA, + *axisB, + *colsize + ); + + processinfo_update_output_stream(processinfo, outimg.ID); + } + INSERT_STD_PROCINFO_COMPUTEFUNC_END + + DEBUG_TRACE_FEXIT(); + return RETURN_SUCCESS; +} + + + +INSERT_STD_FPSCLIfunctions + + + +// Register function in CLI +errno_t +CLIADDCMD_COREMOD_arith__image_unfold() +{ + //CLIcmddata.FPS_customCONFsetup = customCONFsetup; + //CLIcmddata.FPS_customCONFcheck = customCONFcheck; + + INSERT_STD_CLIREGISTERFUNC + + return RETURN_SUCCESS; +} diff --git a/src/COREMOD_arith/image_unfold.h b/src/COREMOD_arith/image_unfold.h new file mode 100644 index 00000000..829d1cfc --- /dev/null +++ b/src/COREMOD_arith/image_unfold.h @@ -0,0 +1,6 @@ +#ifndef COREMOD_ARITH_IMAGE_UNFOLD_H +#define COREMOD_ARITH_IMAGE_UNFOLD_H + +errno_t CLIADDCMD_COREMOD_arith__image_unfold(); + +#endif