TestPatterns.cpp

Go to the documentation of this file.

/* ImLib3D
 * Copyright (c) 2001, ULP-IPB Strasbourg.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
#include<ImLib3D/TestPatterns.hpp>
#include<ImLib3D/MorphologicalOperators.hpp>



void
IP3D::Explode(const Size3D &size,const Vect3Df &C,Field3Df &ima)
{
    ima.Resize(size);
    Field3Df::iteratorXYZ p;
    for(p=ima.begin();p!=ima.end();++p)
    {
        *p=Vect3Df(p.x,p.y,p.z)-C;
    }
}

void
IP3D::Sphere(const Size3D &size,const Vect3Df &C,float r,Image3Df &ima)
{
    ima.Resize(size);
    Image3Df::iteratorXYZ p;
    for(p=ima.begin();p!=ima.end();++p)
    {
        if(C.Dist2(Vect3Df(p.x,p.y,p.z))<r*r){*p=1;}
        else{*p=0;}
    }
}

namespace IP3D
{
namespace ConeHelpers
{
void
NoAA(Image3Df &ima,const Vect3Df &P0,const Vect3Df &P1,float r)
{
    Image3Df::iteratorXYZ p;
    Vect3Df P0P1=P1-P0;
    Vect3Df Uz=P0P1;Uz.Unit();
    float l=P0P1.Norme();
    for(p=ima.begin();p!=ima.end();++p)
    {
        Vect3Df P=Vect3Df(p.x,p.y,p.z);
        // cyl coords:
        float  cz=(P-P0)*Uz;
        Vect3Df N=(P-P0)-cz*Uz;
        float cr=N.Norme();
        if(cz<l && cz>0 && cr<(r*(1-cz/l)) ){*p=1;}
        else{*p=0;}
    }
}
}
}
void
IP3D::Cone(const Size3D &size,const Vect3Df &P0,const Vect3Df &P1,float r,Image3Df &ima,float ar)
{
    ima.Resize(size);
    if(ar==0){ConeHelpers::NoAA(ima,P0,P1,r);return;}
    Image3Df::iteratorXYZ p;
    Vect3Df P0P1=P1-P0;
    Vect3Df Uz=P0P1;Uz.Unit();
    float l=P0P1.Norme();
    float M=M_PI/2.0;
    for(p=ima.begin();p!=ima.end();++p)
    {
        Vect3Df P=Vect3Df(p.x,p.y,p.z);
        // cyl coords:
        float  cz=(P-P0)*Uz;
        if(cz>l){*p=0;continue;}

        Vect3Df N=(P-P0)-cz*Uz;
        float cr=N.Norme();
        float r1=r*(1-cz/l);
        float dr=-(cr-r1)/ar;
        float dz= (cz-0 )/ar;
        bool isBorder=fabs(dz)<1 || fabs(dr)<1;

        if(!isBorder)
        {
            if(cz<l && cz>0 && cr<(r*(1-cz/l)) ){*p=1;}
            else{*p=0;}
        }
        else
        {
            float fz=(dz>0 ? 1 : 0);
            if(fabs(dz)<=1){fz=.5+.5*sin(M*dz);}
            float fr=(dr>0 ? 1 : 0);
            if(fabs(dr)<=1){fr=.5+.5*sin(M*dr);}
            *p=fz*fr;
        }
    }
}
void
IP3D::NoiseUniform(const Size3D &size,Image3Df &ima,float v0,float v1)
{
    ima.Resize(size);
    Image3Df::iteratorXYZ p;
    for(p=ima.begin();p!=ima.end();++p)
    {
        *p=RandomDouble(v0,v1);
    }
}

void
IP3D::NoiseGaussian(const Size3D &size,Image3Df &ima,float average,float variance)
{
    ima.Resize(size);
    Image3Df::iteratorXYZ p;
    double et=sqrt(variance);
    for(p=ima.begin();p!=ima.end();++p)
    {
        double randVal=RandomDouble();
        double r=sqrt(-2 * log(randVal) );
        double angle=RandomDouble(0,2*M_PI);
        *p=average+et*r*cos(angle);
    }
}

void
IP3D::Target(const Size3D &size,Image3Df &ima)
{
    ima.Resize(size);
    Image3Df::iteratorXYZ p;
    for(p=ima.begin();p!=ima.end();++p)
    {
        Vect3Df P(p.x/(float)ima.Width(),p.y/(float)ima.Height(),p.z/(float)ima.Depth()); 
        double d=(P-Vect3Df(.5,.5,.5)).Norme();
        d=d/.5;
        float cx=.8;
        float a=((d-cx)/(1-cx))*M_PI/2.0;
        float f=1;
        if(d>cx && d<1){f=cos(a);f=f*f;}
        if(d>1){f=0;}
        *p=f*sin((.7*128.0)*d*d);
    }
}

void
IP3D::Bump(const Size3D &size,const Vect3Df &center,float height,float gsize,Image3Df &ima)
{
    ima.Resize(size);
    Image3Df::iteratorXYZ p;
    for(p=ima.begin();p!=ima.end();++p)
    {
        Vect3Df pos=p.Pos();
        double r2=pos.Dist2(center);
        *p=height*exp(-r2/(2*gsize*gsize));
    }
}

void
IP3D::Ramp(const Size3D &size,int direction,float v0,float v1,Image3Df &res,int pos0,int pos1)
{
    res.Resize(size);

    if(direction<0 || direction>2){ThrowError("IP3D::Ramp bad direction 0=X 1=Y 2=Z");}
    if(pos0<0){pos0=0;}
    if(pos1<0){pos1=res.SizeV()(direction)-1;}
    for(Image3Df::iteratorXYZ p=res.begin();p!=res.end();++p)
    {
        int pos=p.Pos()(direction);
        if(pos<pos0){*p=v0;}
        else
        if(pos>pos1){*p=v1;}
        else
        {
            *p=v0+(v1-v0)*((pos-pos0)/(float)(pos1-pos0));
        }
    }
}

void
IP3D::Parallellogram(const Size3D &size, const RectZone3Di & rect,Image3Df &ima)
{
    ima.Resize(size);
    Image3Df::iteratorXYZ p;
    
    for (p=ima.begin(); p!=ima.end(); p++)
    {
        if (rect.IsInside(p.x, p.y, p.z)){*p = 1;}
        else{*p = 0;}
    }
}

void
IP3D::RectangularGrid(const Size3D &size, const Size3D& masksize,Image3Df &ima, const Vect3Di* pP0)
{
    ima.Resize(size);
    Vect3Di defP0(0,0,0);
    const Vect3Di &P0=(pP0 ? *pP0 : defP0);
    //Find first intersection
    Vect3Di P = Vect3Di(P0.x%masksize.width, P0.y%masksize.height, P0.z%masksize.depth);

    int x, y, z; //run indexes
    for (x=0;x<ima.Width();x++)
    {
        for (y=0;y<ima.Height();y++)
        {
            for (z=0;z<ima.Depth();z++)
            {
                if (!((P.z+z)%masksize.depth || (P.y+y)%masksize.height)) //on grid z and y coords, we set all x's 1
                    ima(x, y, z)=1;
                else if(!((P.x+x)%masksize.width || (P.y+y)%masksize.height)) //idem for x and y
                    ima(x, y, z)=1;
                else if(!((P.x+x)%masksize.width || (P.z+z)%masksize.depth)) //and x - z
                    ima(x, y, z)=1;
                else
                    ima(x, y, z)=0;
            }
        }
    }
}

void
IP3D::ColoredGrid(const Size3D &size, const Size3D& masksize, Mask3D &ima,StructureElementType structure)
{
    ima.Resize(size);
    // colored by hand
    int x,y,z;
    if (structure == MORPHO_Cross7)
    {
        int lastx = 2;
        int xcount = 0;
        for (x=0;x<ima.Depth();x++, xcount++)
        {
            int lasty = lastx;
            int ycount = 0;
            for (y=0;y<ima.Height();y++, ycount++)
            {
                int lastz = lasty;
                int zcount = 0;
                for (z=0;z<ima.Width();z++, zcount++)
                {
                    if (zcount==masksize.width)
                    {
                        lastz = ima(z,y,x) = (byte) ((lastz==1) ? 2 : 1);
                        zcount = 0;
                    }
                    else
                        ima(z,y,x)=lastz;
                }
                if (ycount == masksize.height)
                {
                    lasty = (lasty==1) ? 2 : 1;
                    ycount = 0;
                }
            }
            if (xcount == masksize.depth)
            {
                lastx = (lastx==1) ? 2 : 1;
                xcount = 0;
            }
        }
    }
    else if (structure == MORPHO_Planar19)
    {
        int actx = 1;
        int xcount = 0;
        for (x=0;x<ima.Depth();x++, xcount++)
        {
            int acty = actx;
            int ycount = 0;
            for (y=0;y<ima.Height();y++, ycount++)
            {
                int actz=acty;
                int zcount = 0;
                for (z=0;z<ima.Width();z++, zcount++)
                {
                    if (zcount==masksize.width)
                    {
                        if (actz==3)
                            actz=4;
                        else if (actz==4)
                            actz=3;
                        else if (actz==1)
                            actz=2;
                        else if (actz==2)
                            actz=1;
                        ima(z, y, x) = (byte)actz;
                        zcount = 0;
                    }
                    else
                        ima(z,y,x)=(byte)actz;
                }
                if (ycount == masksize.height)
                {
                    acty+=2;
                    acty%=4;
                    if (acty==0) acty=4;
                    ycount = 0;
                }
            }
            if (xcount == masksize.depth)
            {
                actx = (actx==1) ? 4 : 1;
                xcount = 0;
            }
        }
    }
    else if (structure == MORPHO_Cubic27)
    {
        int actx = 1;
        int xcount = 0;
        for (x=0;x<ima.Depth();x++, xcount++)
        {
            int acty = actx;
            int ycount = 0;
            for (y=0;y<ima.Height();y++, ycount++)
            {
                int actz=acty;
                int zcount = 0;
                for (z=0;z<ima.Width();z++, zcount++)
                {
                    if (zcount==masksize.width)
                    {
                        if (actz==1||actz==3)
                            actz=(actz==1)?3:1;
                        if (actz==2||actz==6)
                            actz=(actz==2)?6:2;
                        if (actz==4||actz==5)
                            actz=(actz==4)?5:4;
                        if (actz==7||actz==8)
                            actz=(actz==7)?8:7;
                        ima(z, y, x) = (byte)actz;
                        zcount = 0;
                    }
                    else
                        ima(z,y,x)=(byte)actz;
                }
                if (ycount == masksize.height)
                {
                    if (acty==1 || acty==4)
                        acty=(acty==1)?4:1;
                    else if (acty==2 || acty==7)
                        acty=(acty==2)?7:2;
                    ycount = 0;
                }
            }
            if (xcount == masksize.depth)
            {
                actx = (actx==1) ? 2 : 1;
                xcount = 0;
            }
        }
    }
    else 
    {
        ThrowError("TestPatterns::ColoredGrid ERROR No Colored grid for such a structure element implemented!");
    }
}

---