image_find_center.cpp File Reference

#include "core/metadata_vec.h"
#include "core/xmipp_program.h"
#include "core/xmipp_image.h"
#include "core/metadata_sql_operations.h"

Include dependency graph for image_find_center.cpp:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
class	ProgFindCenter2D

Typedefs
typedef unsigned char	BYTE
typedef unsigned short	UWORD
typedef unsigned long	ULONG

Functions
void **	imalloc (int row, int col, int format)
void	imfree (char **image, int row, int col, int format)
float	conv1x (double y, double x)
void	ergrot (double xc0, double yc0, float *zzz)
void	suaviza ()
void	busca ()

Variables
constexpr int	INTFMT = 2
constexpr int	LONGFMT = 3
constexpr int	FLOATFMT = 4
constexpr int	FOURIER = 5
constexpr int	SPIDER = 6
constexpr int	DEF_IT = 5
constexpr int	DEF_DEL = 2
constexpr int	DEF_IN = 2
constexpr int	ALLOC_SIZE = 65000
float	coseno [1281]
BYTE **	imagen
float	r1
float	r2
float	r3
float	cxp1
float	cxp2
float	cxp3
float	cxm1
float	cxm2
float	cxm3
float	rh = 1000.
float	xc0
float	yc0
float	del
float	rbajo
float	ralto
float	zzz0
int	largo
int	lancho
int	indmul
int	in
int	ni
int	mu
int	ncic
int	ir
int	m
int	mu1
int	mu4
int	ntot
int	mt
int	idz
int	ncic2
int	ntot4

Typedef Documentation

BYTE

typedef unsigned char BYTE

Definition at line 48 of file image_find_center.cpp.

ULONG

typedef unsigned long ULONG

Definition at line 50 of file image_find_center.cpp.

UWORD

typedef unsigned short UWORD

Definition at line 49 of file image_find_center.cpp.

Function Documentation

busca()

void busca

(

)

Definition at line 399 of file image_find_center.cpp.

{
    static float f[22][22];
    float pi = 3.141592653;
    short in1, mu5, i, mu3, ind, ind1, indice, iflag, j, ix=0, iy=0, k;
    float h, th, costh, sinth, fi, anc, xk0, yk0, hpi, x, y, z, xx, yy;
    float b, g, d1, e1, ext;
    int count = 0;

    suaviza();
    in = XMIPP_MIN(in, 10);
    mu1 = mu - 1;
    m = 2 * mu;
    mu4 = mu * ncic;
    mt = 2 * m;
    ntot = mt * ncic;
    ntot4 = ntot + mu4 + ncic;
    zzz0 = ntot;
    ncic2 = 2 * ncic;
    in1 = in + 1;
    h = 2. * pi / ntot;
    idz = 2 - ntot;
    th = ncic * h;
    costh = cos(th);
    sinth = sin(th);
    b = 2. / (th * th) * (1. + costh * costh - sin(2. * th) / th);
    g = 4. / (th * th) * (sinth / th - costh);
    d1 = 2. * th / 45.;
    e1 = d1 * sinth * 2.;
    d1 *= sin(2. * th);
    mu5 = mu4 - 1;
    for (i = 1; i <= mu5; i++)
    {
        fi = i * h;
        coseno[i] = sin(fi);
    }
    coseno[mu4] = 1.;
    mu3 = 2 * mu4;
    for (i = 1; i <= mu5; i++)
    {
        coseno[mu3-i] = coseno[i];
        coseno[mu3+i] = -coseno[i];
        coseno[ntot-i] = -coseno[i];
        coseno[ntot+i] = coseno[i];
    }
    coseno[mu3] = 0.;
    coseno[mu3+mu4] = -1.;
    coseno[ntot] = 0.;
    coseno[ntot+mu4] = 1.;
    ind = 2 * in + 1;
    ind1 = ind - 1;
    indice = 0;
    iflag = 0;
e9:
    if (indice >= ni)
        goto e10;
    if (iflag == 2)
        goto e22;
    anc = (int)(3. + in * del);
    xk0 = xc0;
    yk0 = yc0;
    rh = XMIPP_MIN(rh, r2);
    rh = XMIPP_MIN(rh, xk0 - anc - 1.);
    rh = XMIPP_MIN(rh, yk0 - anc - 1.);
    rh = XMIPP_MIN(rh, largo - xk0 - anc);
    rh = XMIPP_MIN(rh, lancho - yk0 - anc);
    ir = (int)((rh - r1) / r3 + 1);
    hpi = h / 2. / pi / ir;
    cxp1 = 2. * b * e1 * hpi;
    cxp2 = -2. * g * d1 * hpi;
    cxp3 = 2. * (g * e1 - b * d1) * hpi;
    cxm1 = (b * b + e1 * e1) * hpi;
    cxm2 = (g * g + d1 * d1) * hpi;
    cxm3 = 2. * (b * g - d1 * e1) * hpi;
    /*    if(iflag == 1) goto 21            */
    if (iflag == 2)
        goto e22;
    x = xc0 - in * del;
    for (i = 1; i <= ind; i++) /* do 5 */
    { y = yc0 - in * del;
        for (j = 1; j <= ind; j++) /*do 4 */
        {   ergrot(x, y, &z);
            /*     printf ("%10.3f%10.3f%10.5f%10.3f%10.3f AA\n", x,y,z,del,rh);*/
            printf(".");
            fflush(stdout);
            z = 100. + indmul * z;
            f[i][j] = z;
            y += del;
        }
        x += del;
    }
    // Introduced by Sjors dd 28.9.2004 to avoid infinite loops
    count++;
    if (count > 1000)
        goto s1;
e23:
    ext = -1000000.;
    for (i = 1; i <= ind; i++) /* do 7 */
        for (j = 1; j <= ind; j++) /* do 7 */
            {   if (f[i][j] > ext)
            {
                ix = i;
                iy = j;
                ext = f[i][j];
            }
        }
    xc0 = xc0 + (ix - in - 1) * del;
    yc0 = yc0 + (iy - in - 1) * del;
    if (ix == 1 || ix == ind || iy == 1 || iy == ind)
        goto e8;
    del /= in;
    iflag = 2;
    indice++;
    goto e9;
e8:
    iflag = 1;
    goto e9;
e22:
    f[1][1] = f[ix-1][iy-1];
    f[ind][1] = f[ix+1][iy-1];
    f[1][ind] = f[ix-1][iy+1];
    f[ind][ind] = f[ix+1][iy+1];
    f[1][in1] = f[ix-1][iy];
    f[in1][1] = f[ix][iy-1];
    f[ind][in1] = f[ix+1][iy];
    f[in1][ind] = f[ix][iy+1];
    f[in1][in1] = f[ix][iy];
    x = xc0 - (in - 1) * del;
    for (i = 2; i < ind1; i++) /* do 11 */
    { y = yc0 - (in - 1) * del;
        for (j = 2;j <= ind1; j++) /* do 12 */
            {   if (i == in1 && j == in1)
                goto e13;
            ergrot(x, y, &z);
            /*     printf ("%10.3f%10.3f%10.5f%10.3f BB \n", x,y,z,del);*/
            printf(".");
            fflush(stdout);
            z = 100. + indmul * z;
            f[i][j] = z;
e13:
            y += del;
        }
        x += del;
    }
    x = xc0 - (in - 1) * del;
    y = yc0 - (in - 1) * del;

    for (k = 2; k <= ind1; k++) /* do 16 */
        { if (k == in1)
            goto e17;
        xx = xc0 - in * del;
        ergrot(xx, y, &z);
        /* printf ("%10.3f%10.3f%10.5f%10.3f CC \n", xx,y,z,del);*/
        printf(".");
        fflush(stdout);
        z = 100. + indmul * z;
        f[1][k] = z;
        yy = yc0 - in * del;
        ergrot(x, yy, &z);
        /* printf ("%10.3f%10.3f%10.5f%10.3f DD \n", xx,y,z,del);*/
        printf(".");
        fflush(stdout);
        z = 100. + indmul * z;
        f[k][1] = z;
        xx = xc0 + in * del;
        ergrot(xx, y, &z);
        /* printf ("%10.3f%10.3f%10.5f%10.3f EE \n", xx,y,z,del);*/
        printf(".");
        fflush(stdout);
        z = 100. + indmul * z;
        f[ind][k] = z;
        yy = yc0 + in * del;
        ergrot(x, yy, &z);
        /* printf ("%10.3f%10.3f%10.5f%10.3f FF\n", x,yy,z,del);*/
        printf(".");
        fflush(stdout);
        z = 100. + indmul * z;
        f[k][ind] = z;
e17:
        x += del;
        y += del;
    }
    goto e23;
e10:
    std::cout << "\nOptimal center coordinates: x= " << yc0-1 << " ,y= " << xc0-1 << " " << std::endl;
    return;
s1:
    yc0=xc0=-1;
    printf("\nNot-converged\n");
    return;
}

conv1x()

float conv1x	(	double	y,
		double	x
	)

Definition at line 224 of file image_find_center.cpp.

{   short i, j;   /* Row and column */
    float intfila1, intfila2; /* Partial row interpolations */
    float escala;   /* Scale factor */

    j = (short int) y;     /* Trunc the x, y coordinates to short */
    i = (short int) x;

    escala =  y - j;
    /* 1st row interpolation */
    intfila1 = imagen[i][j] + escala * ((short)imagen[i][j+1] - (short)imagen[i][j]);
    /* 2nd row interpolation */
    intfila2 = imagen[i+1][j] + escala * ((short)imagen[i+1][j+1] - (short)imagen [i+1][j]);
    /* Column interpolation */
    return intfila1 + (x - i)*(intfila2 - intfila1);
}

ergrot()

void ergrot	(	double	xc0,
		double	yc0,
		float *	zzz
	)

Definition at line 246 of file image_find_center.cpp.

{
    //double xc0, yc0;
    //float *zzz; /* It hides global variable, handle with care */

    static float a[266], b[257];
    double za, zb;
    double xp1, xp2, xp3, xm1, xm2, xm3;
    double axp1, axp2, axp3, axm1, axm2, axm3;
    short i7, iz2, iz1, kv, l1, l2, i1, i, j;
    float r, x, y, zz, ai, bi;
    float bj, aj, z;

    r = r1 - r3;
    axp1 = 0.;
    axm1 = 0.;
    axp2 = 0.;
    axm2 = 0.;
    axp3 = 0.;
    axm3 = 0.;
    for (i7 = 1; i7 <= ir; i7++) /* do 31 */
    { r += r3;
        iz2 = ntot - m + mu4;
        iz1 = iz2 - 1;
        for (i = 1; i <= mu; i++) /* do 13 */
        {   iz1 += idz;
            iz2 += idz;
            za = 0.;
            zb = 0.;
            for (kv = 1; kv <= ncic; kv++) /* do 11 */
            { iz1 += m;
                x = xc0 + r * coseno[iz1];
                y = yc0 + r * coseno[iz1-mu4];
                z = conv1x(y, x);
                iz1 += m;
                x = xc0 + r * coseno[iz1];
                y = yc0 + r * coseno[iz1-mu4];
                zz = conv1x(y, x);
                zb += (z - zz);
                iz2 += m;
                x = xc0 + r * coseno[iz2];
                y = yc0 + r * coseno[iz2-mu4];
                z = conv1x(y, x);
                iz2 += m;
                x = xc0 + r * coseno[iz2];
                y = yc0 + r * coseno[iz2-mu4];
                zz = conv1x(y, x);
                za += (z - zz);
            }
            b[i] = zb;
            a[i] = za;
        }
        xp1 = a[mu];
        xp2 = b[mu];
        xm2 = xp2 * xp2;
        xp3 = xp1 * xp2 * coseno[mu4-ncic];
        xm3 = xp3;
        xp2 = xm2 * coseno[mu4-ncic2];
        xp1 = xp1 * xp1;
        xm1 = xp1;
        for (i = 1; i <= mu1; i++) /* do 14 */
        {   l1 = 4 * i * ncic + mu4;
            l2 = mu4;
            ai = a[i];
            bi = b[i];
            xp1 += (ai * ai * coseno[l1]);
            xp2 += (bi * bi * coseno[l1-ncic2]);
            xp3 += (ai * bi * coseno[l1-ncic]);
            xm1 += (ai * ai);
            xm2 += (bi * bi);
            xm3 += (ai * bi * coseno[mu4+ncic]);
            i1 = i + 1;
            ai = a[i];
            bi = b[i];
            for (j = i1; j <= mu; j++) /* do 15 */
            { l1 += ncic2;
                l2 += ncic2;
                aj = a[j];
                bj = b[j];
                double ajai2=2.0*aj*ai;
                double bjbi2=2.0*bj*bi;
                double aibj=ai*bj;
                double ajbi=aj*bi;
                xp1 += (ajai2 * coseno[l1]);
                xm1 += (ajai2 * coseno[l2]);
                xp2 += (bjbi2 * coseno[l1-ncic2]);
                xm2 += (bjbi2 * coseno[l2]);
                xp3 += ((aibj + ajbi) * coseno[l1-ncic]);
                xm3 += (aibj * coseno[l2-ncic] + ajbi * coseno[l2+ncic]);
            }
        }
        axp1 += xp1 * r;
        axm1 += xm1 * r;
        axp2 += xp2 * r;
        axm2 += xm2 * r;
        axp3 += xp3 * r;
        axm3 += xm3 * r;
    }
    (*zzz) = axp1 * cxp1 + axp2 * cxp2 + axp3 * cxp3 + axm1 * cxm1 + axm2 * cxm2 + axm3 * cxm3;
    (*zzz) /= (zzz0 * ir);
}

imalloc()

void** imalloc	(	int	row,
		int	col,
		int	format
	)

Definition at line 73 of file image_find_center.cpp.

{
    size_t i, j, k;                 /* Counters                             */
    unsigned element_size;          /* Size of element to allocate          */
    unsigned pointer_size;          /* Id. of pointers                      */
    char **temp;                    /* Temporal value to work with          */
    unsigned no_blocks;             /* No. of ALLOC_SIZE blocks to allocate */
    long tot_size;                  /* Total allocation size                */
    unsigned row_alloc;             /* No. of rows to alloc at the same time*/
    unsigned all_size;              /* No. of bits to alloc at one time     */
    unsigned rest_size;             /* Rest of bytes to alloc               */
    unsigned rest_rows;             /* Rest of rows to alloc.               */
    char *aux;                      /* Aux. pointer                         */


    /******************* Assign appropriate value to size flag ******************/

    if (format == FOURIER)
        col += 2;      /* Special treatment for FFT format (see foutrans.c) */

    switch (format)
    {
    case NATURAL:
        element_size = sizeof(BYTE);
        pointer_size = sizeof(BYTE *);
        break;
    case INTFMT:
        element_size = sizeof(UWORD);
        pointer_size = sizeof(UWORD *);
        break;
    case LONGFMT:
        element_size = sizeof(ULONG);
        pointer_size = sizeof(ULONG *);
        break;
    case FLOATFMT:
    case FOURIER:
        element_size = sizeof(float);
        pointer_size = sizeof(float *);
        break;
    default:
        return NULL;
    }

    row_alloc = ALLOC_SIZE / (col * element_size);  /* No. of rows to alloc */
    all_size = row_alloc * col * element_size;
    tot_size = ((long) element_size) * ((long) row) * ((long) col);
    no_blocks = tot_size / all_size;
    rest_size = tot_size - ((long) no_blocks) * ((long) all_size);
    rest_rows = rest_size / (col * element_size);

    /********************* Allocate base pointer ********************************/

    if ((temp = (char **) malloc(row * pointer_size)) == NULL)
        return NULL;               /* Not even this little bit of memory */

    /*********************** Allocate most blocks *******************************/

    j = 0;
    for (i = 0; i < no_blocks; i++)
    {
        if ((aux = (char *)malloc((long)all_size)) == NULL)
        {
            for (j = 0; j < i; j++)
                free(temp[j*row_alloc]);
            free((char *) temp);
            return NULL;
        }
        for (k = 0; k < row_alloc; k++, j++)
            temp [j] = aux + k * col * element_size;
    }

    /*************************** Alloc the last block ***************************/

    if (rest_size != 0)
    {
        if ((aux = (char *)malloc((long)rest_size)) == NULL)
        {
            for (j = 0; j < no_blocks; j++)
                free(temp[j*row_alloc]);
            free(temp);
            return NULL;
        }
        for (k = 0; k < rest_rows; k++, j++)
            temp [j] = aux + k * col * element_size;
    }

    /************************* return OK pointer value  **********************/

    return (void **)temp;
}

imfree()

void imfree	(	char **	image,
		int	row,
		int	col,
		int	format
	)

Definition at line 170 of file image_find_center.cpp.

{
    size_t i;                       /* Counters                             */
    unsigned element_size;          /* Size of element to allocate          */
    unsigned no_blocks;             /* No. of ALLOC_SIZE blocks to allocate */
    long tot_size;                  /* Total allocation size                */
    unsigned row_alloc;             /* No. of rows to alloc at the same time*/
    unsigned all_size;              /* No. of bits to alloc at one time     */

    if (image == NULL)                  /* No allocation at the moment */
        return;

    if (format == FOURIER)
        col += 2;      /* Special treatment for FFT format (see foutrans.c) */

    switch (format)
    {
    case NATURAL:
        element_size = sizeof(BYTE);
        break;
    case INTFMT:
        element_size = sizeof(UWORD);
        break;
    case LONGFMT:
        element_size = sizeof(ULONG);
        break;
    case FLOATFMT:
    case FOURIER:
        element_size = sizeof(float);
        break;
    default:
        return;
    }

    row_alloc = ALLOC_SIZE / (col * element_size);  /* No. of rows to free  */
    all_size = row_alloc * col * element_size;
    tot_size = ((long) element_size) * ((long) row) * ((long) col);
    no_blocks = tot_size / all_size;

    /*************************** Free most blocks *******************************/

    for (i = 0; i < no_blocks; i++)
        if (image [i*row_alloc] != NULL)
            free(image [i*row_alloc]);

    /*************************** Free the last block ****************************/

    if (image [i*row_alloc] != NULL)
        free(image [i*row_alloc]);

    free(image);
}

suaviza()

void suaviza

(

)

Definition at line 348 of file image_find_center.cpp.

{
    unsigned char pix;
    short i, j, k, n;
    long isuma;
    float racua, rbcua, dr, r;

    isuma = 0;
    n = 0;
    i = 0;
    racua = ralto * ralto;
    rbcua = rbajo * rbajo;
    dr = 3.141592653 / (ralto - rbajo);
    for (k = 1; k <= lancho; k++)
    {
        if (k % 16 == 0)
            for (j = 1; j <= largo; j++)
            {
                r = (xc0 - k) * (xc0 - k) + (yc0 - j) * (yc0 - j);
                if (r > racua)
                {
                    isuma += imagen [k][j];
                    n++;
                }
            }
    }
    if (n != 0)
        i = (short int)(((float) isuma) / n + .5);
    m = i;
    pix = i;
    for (k = 1; k <= lancho; k++)
    {
        if (k % 16 == 0)
            for (j = 1; j <= largo; j++)
            {
                r = (xc0 - k) * (xc0 - k) + (yc0 - j) * (yc0 - j);
                if (r < rbcua)
                    continue;
                if (r > racua)
                {
                    imagen[k][j] = pix;
                    continue;
                }
                r = sqrt(r);
                r -= rbajo;
                imagen[k][j] = (unsigned char)((0.5 + 0.5 * cos(dr * r)) * (imagen[k][j] - m) + m);
            }
    }
    printf("\nImage smoothed. Outer mean = %d\n", m);
}

Variable Documentation

ALLOC_SIZE

constexpr int ALLOC_SIZE = 65000

Definition at line 47 of file image_find_center.cpp.

coseno

float coseno[1281]

Definition at line 52 of file image_find_center.cpp.

cxm1

float cxm1

Definition at line 54 of file image_find_center.cpp.

cxm2

float cxm2

Definition at line 54 of file image_find_center.cpp.

cxm3

float cxm3

Definition at line 54 of file image_find_center.cpp.

cxp1

float cxp1

Definition at line 54 of file image_find_center.cpp.

cxp2

float cxp2

Definition at line 54 of file image_find_center.cpp.

cxp3

float cxp3

Definition at line 54 of file image_find_center.cpp.

DEF_DEL

constexpr int DEF_DEL = 2

Definition at line 45 of file image_find_center.cpp.

DEF_IN

constexpr int DEF_IN = 2

Definition at line 46 of file image_find_center.cpp.

DEF_IT

constexpr int DEF_IT = 5

Definition at line 44 of file image_find_center.cpp.

del

float del

Definition at line 54 of file image_find_center.cpp.

FLOATFMT

constexpr int FLOATFMT = 4

Definition at line 41 of file image_find_center.cpp.

FOURIER

constexpr int FOURIER = 5

Definition at line 42 of file image_find_center.cpp.

idz

int idz

Definition at line 57 of file image_find_center.cpp.

imagen

BYTE** imagen

Definition at line 53 of file image_find_center.cpp.

in

int in

Definition at line 56 of file image_find_center.cpp.

indmul

int indmul

Definition at line 56 of file image_find_center.cpp.

INTFMT

constexpr int INTFMT = 2

Definition at line 39 of file image_find_center.cpp.

ir

int ir

Definition at line 57 of file image_find_center.cpp.

lancho

int lancho

Definition at line 56 of file image_find_center.cpp.

largo

int largo

Definition at line 56 of file image_find_center.cpp.

LONGFMT

constexpr int LONGFMT = 3

Definition at line 40 of file image_find_center.cpp.

m

if m

Initial value:

{
        goto L45

Definition at line 57 of file image_find_center.cpp.

mt

int mt

Definition at line 57 of file image_find_center.cpp.

mu

int mu

Definition at line 56 of file image_find_center.cpp.

mu1

int mu1

Definition at line 57 of file image_find_center.cpp.

mu4

int mu4

Definition at line 57 of file image_find_center.cpp.

ncic

int ncic

Definition at line 56 of file image_find_center.cpp.

ncic2

int ncic2

Definition at line 57 of file image_find_center.cpp.

ni

int ni

Definition at line 56 of file image_find_center.cpp.

ntot

int ntot

Definition at line 57 of file image_find_center.cpp.

ntot4

int ntot4

Definition at line 57 of file image_find_center.cpp.

r1

float r1

Definition at line 54 of file image_find_center.cpp.

r2

float r2

Definition at line 54 of file image_find_center.cpp.

r3

float r3

Definition at line 54 of file image_find_center.cpp.

ralto

float ralto

Definition at line 54 of file image_find_center.cpp.

rbajo

float rbajo

Definition at line 54 of file image_find_center.cpp.

rh

float rh = 1000.

Definition at line 54 of file image_find_center.cpp.

SPIDER

constexpr int SPIDER = 6

Definition at line 43 of file image_find_center.cpp.

xc0

float xc0

Definition at line 54 of file image_find_center.cpp.

yc0

float yc0

Definition at line 54 of file image_find_center.cpp.

zzz0

float zzz0

Definition at line 54 of file image_find_center.cpp.