diff --git a/build_linux64/libamscppimglib4.linux64.a b/build_linux64/libamscppimglib4.linux64.a
index d683232..3605aef 100644
Binary files a/build_linux64/libamscppimglib4.linux64.a and b/build_linux64/libamscppimglib4.linux64.a differ
diff --git a/build_linux64/objstore/amscppimglib4_amsfloatimage.o b/build_linux64/objstore/amscppimglib4_amsfloatimage.o
new file mode 100644
index 0000000..dd9002e
Binary files /dev/null and b/build_linux64/objstore/amscppimglib4_amsfloatimage.o differ
diff --git a/build_linux64/objstore/amscppimglib4_intlutil.o b/build_linux64/objstore/amscppimglib4_intlutil.o
index 66c10c5..f2368ff 100644
Binary files a/build_linux64/objstore/amscppimglib4_intlutil.o and b/build_linux64/objstore/amscppimglib4_intlutil.o differ
diff --git a/include/amscppimglib4/amscppimglib4.hpp b/include/amscppimglib4/amscppimglib4.hpp
index 5dcd49e..fd2621a 100644
--- a/include/amscppimglib4/amscppimglib4.hpp
+++ b/include/amscppimglib4/amscppimglib4.hpp
@@ -160,6 +160,53 @@ namespace ams
         amsbitplane rescale(int _Nx, int _Ny); //todo
     };
 
+    class amsfloatimage
+    {
+        public:
+        int Nx,Ny; //image extent
+        //int Nc; //number of color-planes  (1,3,4)
+        //Assume we're talking about an RGBA image internally, greyscale is what bitplane is for
+        float *data; //[(x + width*y)*4]
+
+        int& width;  //aliases to Nx,Ny
+        int& height;
+
+        amsfloatimage();
+        ~amsfloatimage();
+        amsfloatimage(const amsfloatimage& other);
+        amsfloatimage(amsfloatimage&& other) noexcept;
+        amsfloatimage& operator=(const amsfloatimage& other);
+        amsfloatimage& operator=(amsfloatimage&& other) noexcept;
+
+        int resize(int _Nx, int _Ny);
+
+        amsfloatimage subimage(int I0, int J0, int I1, int J1) const;
+        amsfloatimage transpose() const;
+        amsfloatimage rotcw() const;
+        amsfloatimage rotccw() const;
+        amsfloatimage flipx() const;
+        amsfloatimage flipy() const;
+
+        amsfloatpixel get_pixel(int I,int J) const;
+        int set_pixel(int I,int J,const amsfloatpixel pix);
+        int set_pixel(int I,int J,float R, float G, float B, float A);
+        float& operator[](int ind);
+        const float& operator[](int ind) const;
+        float& operator()(int Nc, int I, int J);
+        const float& operator()(int Nc, int I, int J) const;
+
+        int apply_image(int x0, int y0, const amsfloatimage *img);
+
+        void clear();
+        void setall(amsfloatpixel color);
+
+        //don't implement this yet
+        amsfloatpixel interpolate(float x, float y) const;
+
+        // //rescales the image with linear interpolation
+        amsfloatimage rescale(int nnx, int nny);
+    };
+
 
 
 
diff --git a/include/amscppimglib4/amscppimglib4_intlutil.hpp b/include/amscppimglib4/amscppimglib4_intlutil.hpp
index 5f48aa2..49642a8 100644
--- a/include/amscppimglib4/amscppimglib4_intlutil.hpp
+++ b/include/amscppimglib4/amscppimglib4_intlutil.hpp
@@ -109,6 +109,50 @@ float mod(float x, float n);
 int32_t mod(int32_t x, int32_t n);
 int64_t mod(int64_t x, int64_t n);
 
+
+void amsfloatimage_region_set(
+    float *data,
+    int Nx, int Ny,
+    int x0, int y0,
+    int x1, int y1,
+    amsfloatpixel val
+);
+
+void amsfloatimage_region_copy(
+    float *datato,
+    int Nxto,
+    int Nyto,
+    const float *datafrom,
+    int Nxfrom,
+    int Nyfrom,
+    int offsetx,
+    int offsety
+);
+
+void amsfloatimage_region_castcopy(
+    uint8_t *datato,
+    int Nxto,
+    int Nyto,
+    const float *datafrom,
+    int Nxfrom,
+    int Nyfrom,
+    int offsetx,
+    int offsety
+);
+
+void amsfloatimage_region_castcopy(
+    float *datato,
+    int Nxto,
+    int Nyto,
+    const uint8_t *datafrom,
+    int Nxfrom,
+    int Nyfrom,
+    int offsetx,
+    int offsety
+);
+
+
+
 }; //end namespace imglib4
 }; //end namespace ams
 
diff --git a/src/amscppimglib4/amscppimglib4_amsfloatimage.cpp b/src/amscppimglib4/amscppimglib4_amsfloatimage.cpp
new file mode 100644
index 0000000..b640ac2
--- /dev/null
+++ b/src/amscppimglib4/amscppimglib4_amsfloatimage.cpp
@@ -0,0 +1,769 @@
+#include <amscppimglib4/amscppimglib4.hpp>
+#include <amscppimglib4/amscppimglib4_intlutil.hpp>
+#include <amscimglib4/amscimglib4.h>
+
+namespace ams
+{
+
+    amsfloatimage::amsfloatimage() : Nx(0), Ny(0),data(NULL),width(Nx),height(Ny)
+    {
+        return;
+    }
+
+    amsfloatimage::~amsfloatimage()
+    {
+        Nx = 0;
+        Ny = 0;
+        if(data!=NULL) {delete[] data; data=NULL;}
+        return;
+    }
+
+    int amsfloatimage::resize(int _Nx, int _Ny)
+    {
+        int ret = amsimage_success;
+        
+        float *newdata = NULL;
+
+        _Nx = (_Nx<0) ? 0 : _Nx;
+        _Ny = (_Ny<0) ? 0 : _Ny;
+
+        if(_Nx == Nx && _Ny == Ny)
+        {
+            return ret; //no resize necessary
+        }
+
+        if(_Nx==0 || _Ny == 0)
+        {
+            //zero size image
+            if(data!=NULL) {delete[] data; data=NULL;}
+            Nx = 0;
+            Ny = 0;
+            return ret;
+        }
+        
+        newdata = new(std::nothrow) float[4*_Nx*_Ny];
+        if(newdata==NULL)
+        {
+            ret = amsimage_failure;
+            return ret;
+        }
+
+        imglib4::amsfloatimage_region_set(
+            newdata,
+            _Nx,_Ny,
+            0,0,_Nx,_Ny,
+            amsfloatpixel(0,0,0,0)
+        );
+
+        if(data!=NULL)
+        {
+            imglib4::amsfloatimage_region_copy(
+                newdata,
+                _Nx,_Ny,
+                data,
+                Nx,Ny,
+                0,0
+            );
+        }
+
+        if(data!=NULL) {delete[] data; data=NULL;}
+        data = newdata;
+        Nx = _Nx;
+        Ny = _Ny;
+
+        return ret;
+    }
+
+    amsfloatimage::amsfloatimage(const amsfloatimage& other) : 
+        Nx(0), Ny(0),data(NULL),width(Nx),height(Ny)
+    {
+        int res;
+        // Nx = 0;
+        // Ny = 0;
+        // data = NULL;
+
+        if(this!=&other)
+        {
+            res = this->resize(other.Nx,other.Ny);
+            if(res==amsimage_success)
+            {
+                imglib4::amsfloatimage_region_copy(
+                    data, Nx, Ny,
+                    other.data, other.Nx, other.Ny,
+                    0,0
+                );
+            }
+        }
+    }
+
+    amsfloatimage::amsfloatimage(amsfloatimage&& other) noexcept :
+        Nx(0), Ny(0),data(NULL),width(Nx),height(Ny)
+    {
+        int res;
+        // Nx = 0;
+        // Ny = 0;
+        // data = NULL;
+
+        if(this!=&other)
+        {
+            this->Nx = other.Nx;
+            this->Ny = other.Ny;
+            this->data = other.data;
+
+            other.Nx = 0;
+            other.Ny = 0;
+            other.data = NULL;
+        }
+
+        return;
+    }
+
+    amsfloatimage& amsfloatimage::operator=(const amsfloatimage& other)
+    {
+        int res;
+
+        if(this!=&other)
+        {
+            res = this->resize(other.Nx,other.Ny);
+            if(res==amsimage_success)
+            {
+                imglib4::amsfloatimage_region_copy(
+                    data, Nx, Ny,
+                    other.data, other.Nx, other.Ny,
+                    0,0
+                );
+            }
+        }
+
+        return *this;
+    }
+    
+    amsfloatimage& amsfloatimage::operator=(amsfloatimage&& other) noexcept
+    {
+        if(this!=&other)
+        {
+            if(this->data!=NULL) {delete[] this->data; this->data=NULL;}
+            this->Nx = 0;
+            this->Ny = 0;
+
+            this->Nx = other.Nx;
+            this->Ny = other.Ny;
+            this->data = other.data;
+
+            other.Nx = 0;
+            other.Ny = 0;
+            other.data = NULL;
+        }
+
+        return *this;
+    }
+
+    amsfloatimage amsfloatimage::subimage(int I0, int J0, int I1, int J1) const
+    {
+        amsfloatimage ret;
+        int _Nx,_Ny;
+        _Nx = ((I1-I0) < 0) ? 0 : I1-I0;
+        _Ny = ((J1-J0) < 0) ? 0 : J1-J0;
+
+        ret.resize(_Nx,_Ny);
+
+        imglib4::amsfloatimage_region_copy(
+            ret.data, ret.Nx, ret.Ny,
+            this->data, Nx,Ny,
+            I0,J0
+        );
+
+        return ret;
+    }
+
+    void amsfloatimage_transpose_tf(
+        int threadnum,
+        int nthreads,
+        float *datato,
+        const float *datafrom,
+        int Nx, int Ny
+    )
+    {
+        int64_t I,I0,I1,Is,N,Ia,Ib,Ix,Iy;
+        N = Nx*Ny;
+        Is = N/nthreads;
+        I0 = (threadnum)*Is;
+        I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+        for(I=I0;I<I1;I++)
+        {
+            Ix = I%Nx;
+            Iy = I/Nx;
+            Ia = Iy + Ny*Ix;
+            Ib = Ix+Nx*Iy;
+
+            datato[0 + 4*Ia] = datafrom[0 + 4*Ib];
+            datato[1 + 4*Ia] = datafrom[1 + 4*Ib];
+            datato[2 + 4*Ia] = datafrom[2 + 4*Ib];
+            datato[3 + 4*Ia] = datafrom[3 + 4*Ib];
+        }
+        return;
+    }
+
+    amsfloatimage amsfloatimage::transpose() const
+    {
+        amsfloatimage ret;
+        int64_t N;
+
+        ret.resize(Ny,Nx);
+        N = Nx*Ny;
+
+        imglib4::threaded_execute(
+            amsfloatimage_transpose_tf,
+            N,
+            ret.data,this->data,Nx,Ny
+        );
+
+        return ret;
+    }
+
+    void amsfloatimage_rotcw_tf(
+        int threadnum,
+        int nthreads,
+        float *datato,
+        const float *datafrom,
+        int Nx, int Ny
+    )
+    {
+        int64_t I,I0,I1,Is,N,Ia,Ib,Ix,Iy;
+        N = Nx*Ny;
+        Is = N/nthreads;
+        I0 = (threadnum)*Is;
+        I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+        for(I=I0;I<I1;I++)
+        {
+            Ix = I%Nx;
+            Iy = I/Nx;
+            Ia = Iy + Ny*Ix;
+            Ib = Ix+Nx*(Ny-Iy-1);
+
+            datato[0 + 4*Ia] = datafrom[0 + 4*Ib];
+            datato[1 + 4*Ia] = datafrom[1 + 4*Ib];
+            datato[2 + 4*Ia] = datafrom[2 + 4*Ib];
+            datato[3 + 4*Ia] = datafrom[3 + 4*Ib];
+        }
+        return;
+    }
+
+    amsfloatimage amsfloatimage::rotcw() const
+    {
+        amsfloatimage ret;
+        int64_t N;
+
+        ret.resize(Ny,Nx);
+        N = Nx*Ny;
+
+        imglib4::threaded_execute(
+            amsfloatimage_rotcw_tf,
+            N,
+            ret.data,this->data,Nx,Ny
+        );
+
+        return ret;
+    }
+
+    void amsfloatimage_rotccw_tf(
+        int threadnum,
+        int nthreads,
+        float *datato,
+        const float *datafrom,
+        int Nx, int Ny
+    )
+    {
+        int64_t I,I0,I1,Is,N,Ia,Ib,Ix,Iy;
+        N = Nx*Ny;
+        Is = N/nthreads;
+        I0 = (threadnum)*Is;
+        I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+        for(I=I0;I<I1;I++)
+        {
+            Ix = I%Nx;
+            Iy = I/Nx;
+            Ia = Iy + Ny*Ix;
+            Ib = (Nx-Ix-1)+Nx*Iy;
+
+            datato[0 + 4*Ia] = datafrom[0 + 4*Ib];
+            datato[1 + 4*Ia] = datafrom[1 + 4*Ib];
+            datato[2 + 4*Ia] = datafrom[2 + 4*Ib];
+            datato[3 + 4*Ia] = datafrom[3 + 4*Ib];
+        }
+        return;
+    }
+
+    amsfloatimage amsfloatimage::rotccw() const
+    {
+        amsfloatimage ret;
+        int64_t N;
+
+        ret.resize(Ny,Nx);
+        N = Nx*Ny;
+
+        imglib4::threaded_execute(
+            amsfloatimage_rotccw_tf,
+            N,
+            ret.data,this->data,Nx,Ny
+        );
+
+        return ret;
+    }
+
+    void amsfloatimage_flipx_tf(
+        int threadnum,
+        int nthreads,
+        float *datato,
+        const float *datafrom,
+        int Nx, int Ny
+    )
+    {
+        int64_t I,I0,I1,Is,N,Ia,Ib,Ix,Iy;
+        N = Nx*Ny;
+        Is = N/nthreads;
+        I0 = (threadnum)*Is;
+        I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+        for(I=I0;I<I1;I++)
+        {
+            Ix = I%Nx;
+            Iy = I/Nx;
+            Ia = Ix + Nx*Iy;
+            Ib = (Nx-Ix-1)+Nx*Iy;
+
+            datato[0 + 4*Ia] = datafrom[0 + 4*Ib];
+            datato[1 + 4*Ia] = datafrom[1 + 4*Ib];
+            datato[2 + 4*Ia] = datafrom[2 + 4*Ib];
+            datato[3 + 4*Ia] = datafrom[3 + 4*Ib];
+        }
+        return;
+    }
+
+    amsfloatimage amsfloatimage::flipx() const
+    {
+        amsfloatimage ret;
+        int64_t N;
+
+        ret.resize(Nx,Ny);
+        N = Nx*Ny;
+
+        imglib4::threaded_execute(
+            amsfloatimage_flipx_tf,
+            N,
+            ret.data,this->data,Nx,Ny
+        );
+
+        return ret;
+    }
+
+    void amsfloatimage_flipy_tf(
+        int threadnum,
+        int nthreads,
+        float *datato,
+        const float *datafrom,
+        int Nx, int Ny
+    )
+    {
+        int64_t I,I0,I1,Is,N,Ia,Ib,Ix,Iy;
+        N = Nx*Ny;
+        Is = N/nthreads;
+        I0 = (threadnum)*Is;
+        I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+        for(I=I0;I<I1;I++)
+        {
+            Ix = I%Nx;
+            Iy = I/Nx;
+            Ia = Ix + Nx*Iy;
+            Ib = Ix+Nx*(Ny-Iy-1);
+
+            datato[0 + 4*Ia] = datafrom[0 + 4*Ib];
+            datato[1 + 4*Ia] = datafrom[1 + 4*Ib];
+            datato[2 + 4*Ia] = datafrom[2 + 4*Ib];
+            datato[3 + 4*Ia] = datafrom[3 + 4*Ib];
+        }
+        return;
+    }
+
+    amsfloatimage amsfloatimage::flipy() const
+    {
+        amsfloatimage ret;
+        int64_t N;
+
+        ret.resize(Nx,Ny);
+        N = Nx*Ny;
+
+        imglib4::threaded_execute(
+            amsfloatimage_flipy_tf,
+            N,
+            ret.data,this->data,Nx,Ny
+        );
+
+        return ret;
+    }
+
+
+    amsfloatpixel amsfloatimage::get_pixel(int I,int J) const
+    {
+        amsfloatpixel ret;
+        if(I>=0 && I<Nx && J>=0 && J<Ny)
+        {
+            ret.R = data[0 + 4*(I + Nx*J)];
+            ret.G = data[1 + 4*(I + Nx*J)];
+            ret.B = data[2 + 4*(I + Nx*J)];
+            ret.A = data[3 + 4*(I + Nx*J)];
+        }
+        return ret;
+    }
+
+    int amsfloatimage::set_pixel(int I,int J,const amsfloatpixel pix)
+    {
+        int ret = amsimage_failure;
+        if(I>=0 && I<Nx && J>=0 && J<Ny)
+        {
+            ret = amsimage_success;
+            data[0 + 4*(I + Nx*J)] = pix.R;
+            data[1 + 4*(I + Nx*J)] = pix.G;
+            data[2 + 4*(I + Nx*J)] = pix.B;
+            data[3 + 4*(I + Nx*J)] = pix.A;
+        }
+        return ret;
+    }
+
+    int amsfloatimage::set_pixel(int I,int J,float R, float G, float B, float A)
+    {  
+        return this->set_pixel(I,J,amsfloatpixel(R,G,B,A));
+    }
+
+    float& amsfloatimage::operator[](int ind)
+    {
+        return data[ind];
+    }
+
+    const float& amsfloatimage::operator[](int ind) const
+    {
+        return data[ind];
+    }
+
+    float& amsfloatimage::operator()(int Nc, int I, int J)
+    {
+        return data[Nc + 4*(I + Nx*J)];
+    }
+
+    const float& amsfloatimage::operator()(int Nc, int I, int J) const
+    {
+        return data[Nc + 4*(I + Nx*J)];
+    }
+
+
+    int read_image(const char *fname, amsfloatimage* image)
+    {
+        int ret = amsimage_success;
+        int res;
+        amscimglib4_image *im2 = NULL;
+
+        if(image==NULL)
+        {
+            ret = amsimage_failure;
+            printf("read_image: Error: image pointer is null.\n");
+            return ret;
+        }
+
+        res = amscimglib4_image_new(&im2,1,1);
+        if(res!=amscimglib4_success)
+        {
+            ret = amsimage_failure;
+            printf("read_image: Error: c image struct failed to allocate.\n");
+            return ret;
+        }
+
+        amscimglib4_readimage(fname,im2);
+        //you can NOT move buffers. The im2 buffers are created with malloc, not new. No pseudo-move-semantics for you!
+
+        //copy buffers
+        res =  image->resize(im2->sizex,im2->sizey);
+        if(res!=amsimage_success)
+        {
+            ret = amsimage_failure;
+            amscimglib4_image_delete(&im2);
+            printf("read_image: Error: c++ image failed to allocate.\n");
+            return ret;
+        }
+
+        //both structures have the same memory layout, so the internal copy function works
+        imglib4::amsfloatimage_region_castcopy(
+            image->data,image->Nx,image->Ny,
+            im2->data,im2->sizex,im2->sizey,
+            0,0
+        );
+
+
+        amscimglib4_image_delete(&im2);
+
+        return ret;
+    }
+
+    int write_image(const char *fname, amsfloatimage* image)
+    {
+        int ret = amsimage_success;
+        int res;
+        amscimglib4_image *im2 = NULL;
+
+        if(image==NULL)
+        {
+            ret = amsimage_failure;
+            printf("write_image: Error: image pointer is null.\n");
+            return ret;
+        }
+
+        res = amscimglib4_image_new(&im2,image->Nx,image->Ny);
+        if(res!=amscimglib4_success)
+        {
+            ret = amsimage_failure;
+            printf("write_image: Error: c image struct failed to allocate.\n");
+            return ret;
+        }
+
+        //both structures have the same memory layout, so the internal copy function works
+        imglib4::amsfloatimage_region_castcopy(
+            im2->data,im2->sizex,im2->sizey,
+            image->data,image->Nx,image->Ny,
+            0,0
+        );
+
+        amscimglib4_writeimage(fname,im2);
+
+        amscimglib4_image_delete(&im2);
+        return ret;
+    }
+
+    void amsfloatimage::clear()
+    {
+        this->setall(amsfloatpixel(0,0,0,0));
+    }
+
+    void amsfloatimage_setall_tf(
+        int threadnum,
+        int nthreads,
+        amsfloatimage *img,
+        const amsfloatpixel color
+    )
+    {
+        int64_t I,I0,I1,Is,N,Ix,Iy;
+
+        N = img->Nx*img->Ny;
+        Is = N/nthreads; Is = (Is<1) ? 1 : Is;
+        I0 = Is*threadnum;
+        I1 = (threadnum<(nthreads-1)) ? Is*(threadnum+1) : N;
+
+        for(I=I0;I<I1;I++)
+        {
+            Ix = I%img->Nx;
+            Iy = I/img->Nx;
+
+            img->data[0 + 4*(Ix + img->Nx*Iy)] = color.R;
+            img->data[1 + 4*(Ix + img->Nx*Iy)] = color.G;
+            img->data[2 + 4*(Ix + img->Nx*Iy)] = color.B;
+            img->data[3 + 4*(Ix + img->Nx*Iy)] = color.A;
+        }
+    }
+
+    void amsfloatimage::setall(amsfloatpixel color)
+    {
+        imglib4::threaded_execute(
+            amsfloatimage_setall_tf,
+            this->Nx*this->Ny,
+            this,
+            color
+        );
+        return;
+    }
+
+    void amsfloatimage_apply_image_tf(
+        int threadnum,
+        int nthreads,
+        amsfloatimage *imgto,
+        const amsfloatimage *imgfrom,
+        int x0, int y0
+    )
+    {
+        int dx,dy;
+        int64_t N;
+        int64_t I,I0,I1,Is,Ix,Iy,Ia,Ib;
+        double r1,g1,b1,a1;
+        double r2,g2,b2,a2;
+        double r3,g3,b3,a3;
+
+        dx = (imgfrom->Nx < (imgto->Nx-x0)) ? imgfrom->Nx : imgto->Nx-x0;
+        dx = (dx<0) ? 0 : dx;
+        dy = (imgfrom->Ny < (imgto->Ny-y0)) ? imgfrom->Ny : imgto->Ny-y0;
+        dy = (dy<0) ? 0 : dy;
+        N = dx*dy;
+
+        Is = N/nthreads;
+        I0 = threadnum*Is;
+        I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+        for(I=I0;I<I1;I++)
+        {
+            Ix = I%dx;
+            Iy = I/dx;
+
+            Ia = (Ix + x0) + (Iy + y0)*imgto->Ny;
+            Ib = Ix + Iy*imgfrom->Ny;
+
+            r1 = imgto->data[0 + 4*Ia];
+            g1 = imgto->data[1 + 4*Ia];
+            b1 = imgto->data[2 + 4*Ia];
+            a1 = imgto->data[3 + 4*Ia];
+            
+            r2 = imgfrom->data[0 + 4*Ib];
+            g2 = imgfrom->data[1 + 4*Ib];
+            b2 = imgfrom->data[2 + 4*Ib];
+            a2 = imgfrom->data[3 + 4*Ib];
+
+            r3 = r1 + r2*(1.0-a2);
+            g3 = g1 + g2*(1.0-a2);
+            b3 = b1 + b2*(1.0-a2);
+            a3 = 1.0-(1.0-a1)*(1.0-a2);
+
+            imgto->data[0 + 4*Ia] = r3;
+            imgto->data[1 + 4*Ia] = g3;
+            imgto->data[2 + 4*Ia] = b3;
+            imgto->data[3 + 4*Ia] = a3;
+        }
+
+        return;
+    }
+
+    int amsfloatimage::apply_image(int x0, int y0, const amsfloatimage *img)
+    {
+        int ret = amsimage_success;
+
+        int dx,dy;
+        int64_t N;
+
+        if(img==NULL)
+        {
+            ret = amsimage_failure;
+            return ret;
+        }
+
+        dx = (img->Nx < (this->Nx-x0)) ? img->Nx : this->Nx-x0;
+        dx = (dx<0) ? 0 : dx;
+        dy = (img->Ny < (this->Ny-y0)) ? img->Ny : this->Ny-y0;
+        dy = (dy<0) ? 0 : dy;
+        N = dx*dy;
+
+        imglib4::threaded_execute(
+            amsfloatimage_apply_image_tf, N,
+            this,img,x0,y0
+        );
+
+        return ret;
+    }
+
+    amsfloatpixel amsfloatimage::interpolate(float x, float y) const
+    {
+        amsfloatpixel ret = amsfloatpixel(0.0,0.0,0.0,0.0);
+
+        float Nxf,Nyf;
+        //just do simple (0,N-1) interpolation for now
+        int xi,yi;
+        float xif,yif,xr,yr;
+        float w00,w01,w10,w11;
+        amsfloatpixel p00,p01,p10,p11;
+
+        Nxf = (float)Nx;
+        Nyf = (float)Ny;
+
+        if(x>-0.5 && x<Nxf-0.5 && y>-0.5 && y<Nyf-0.5)
+        {
+            xr = imglib4::mod(x,1.0f);
+            yr = imglib4::mod(y,1.0f);
+            xif = x - xr;
+            yif = y - yr;
+            if(xif<0.0f) {xr = 1.0f;}
+            if(xif>=(Nxf-1.0f)) {xr = 0.0f;}
+            if(yif<0.0f) {yr = 1.0;}
+            if(yif>=(Nyf-1.0f)) {yr = 0.0f;}
+            xi = (int)xif;
+            yi = (int)yif;
+            w00 = (1.0-xr)*(1.0-yr);
+            w10 = (xr)*(1.0-yr);
+            w01 = (1.0-xr)*(yr);
+            w11 = (xr)*(yr);
+
+            if(xi<0) xi = 0;
+            if(xi>Nx-2) xi = Nx-2;
+            if(yi<0) yi = 0;
+            if(yi>Ny-2) yi = Ny-2;
+
+            p00 = get_pixel(xi,yi);
+            p10 = get_pixel(xi+1,yi);
+            p01 = get_pixel(xi,yi+1);
+            p11 = get_pixel(xi+1,yi+1);
+
+            ret = p00*w00 + p10*w10 + p01*w01 + p11*w11;
+        }
+
+        return ret;
+    }
+
+    void amsfloatimage_rescale_tf(
+        int threadnum,
+        int nthreads,
+        amsfloatimage *imgto,
+        const amsfloatimage *imgfrom
+    )
+    {
+        int Nx = imgto->Nx;
+        int Ny = imgto->Ny;
+        int64_t I,I0,I1,Is,N,Ix,Iy;
+        amsfloatpixel p;
+
+        int xs,ys;
+
+        N = Nx*Ny;
+        Is = N/nthreads; if(Is<1) Is = 1;
+        I0 = threadnum*Is;
+        I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+        for(I=I0;I<I1;I++)
+        {
+            Ix = I%Nx;
+            Iy = I/Nx;
+
+            xs = (float)Ix*((float)(imgfrom->Nx-1))/((float)(imgto->Nx-1));
+            ys = (float)Iy*((float)(imgfrom->Ny-1))/((float)(imgto->Ny-1));
+
+            p = imgfrom->interpolate(xs,ys);
+            imgto->set_pixel(Ix,Iy,p);
+        }
+    }
+
+    amsfloatimage amsfloatimage::rescale(int nnx, int nny)
+    {
+        amsfloatimage ret;
+        int res;
+
+        res = ret.resize(nnx,nny);
+        if(res!=amsimage_success)
+        {
+            printf("amsfloatimage::rescale: error, could not allocate return image.\n");
+            return ret;
+        }
+
+        imglib4::threaded_execute(
+            amsfloatimage_rescale_tf, nnx*nny,
+            &ret,this
+        );
+
+        return ret;
+    }
+
+};
\ No newline at end of file
diff --git a/src/amscppimglib4/amscppimglib4_intlutil.cpp b/src/amscppimglib4/amscppimglib4_intlutil.cpp
index ffc2751..d0615e8 100644
--- a/src/amscppimglib4/amscppimglib4_intlutil.cpp
+++ b/src/amscppimglib4/amscppimglib4_intlutil.cpp
@@ -357,5 +357,363 @@ int64_t mod(int64_t x, int64_t n)
     return x;
 }
 
+// Floating point region copy/set
+
+void amsfloatimage_region_copy_tf(
+    int threadnum,
+    int nthreads,
+    float *datato,
+    int Nxto,
+    int Nyto,
+    const float *datafrom,
+    int Nxfrom,
+    int Nyfrom,
+    int offsetx,
+    int offsety
+)
+{
+    int64_t I,I0,I1,Is,N,Ix,Iy;
+    int dx,dy;
+
+    dx = Nxfrom;
+    dx = (dx>(Nxto-offsetx)) ? (Nxto-offsetx) : dx;
+    dx = (dx<0) ? 0 : dx;
+
+    dy = Nyfrom;
+    dy = (dy>(Nyto-offsety)) ? (Nyto-offsety) : dy;
+    dy = (dy<0) ? 0 : dy;
+
+    N = dx*dy;
+
+    Is = N/nthreads; Is = (Is<1) ? 1 : N;
+    I0 = (threadnum)*Is;
+    I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+    for(I=I0;I<I1;I++)
+    {
+        Ix = I%dx;
+        Iy = I/dx;
+
+        datato[0 + 4*(Ix+offsetx) + 4*Nxto*(Iy+offsety)] = 
+        datafrom[0 + 4*Ix + 4*Nxfrom*Iy];
+
+        datato[1 + 4*(Ix+offsetx) + 4*Nxto*(Iy+offsety)] = 
+        datafrom[1 + 4*Ix + 4*Nxfrom*Iy];
+        
+        datato[2 + 4*(Ix+offsetx) + 4*Nxto*(Iy+offsety)] = 
+        datafrom[2 + 4*Ix + 4*Nxfrom*Iy];
+        
+        datato[3 + 4*(Ix+offsetx) + 4*Nxto*(Iy+offsety)] = 
+        datafrom[3 + 4*Ix + 4*Nxfrom*Iy];
+    }
+
+    return;
+}
+
+void amsfloatimage_region_copy(
+    float *datato,
+    int Nxto,
+    int Nyto,
+    const float *datafrom,
+    int Nxfrom,
+    int Nyfrom,
+    int offsetx,
+    int offsety
+)
+{
+    int dx,dy;
+    int64_t N;
+
+    if(datato==NULL) return;
+    if(datafrom==NULL) return;
+
+    dx = Nxfrom;
+    dx = (dx>(Nxto-offsetx)) ? (Nxto-offsetx) : dx;
+    dx = (dx<0) ? 0 : dx;
+
+    dy = Nyfrom;
+    dy = (dy>(Nyto-offsety)) ? (Nyto-offsety) : dy;
+    dy = (dy<0) ? 0 : dy;
+
+    N = dx*dy;
+
+    threaded_execute(
+        amsfloatimage_region_copy_tf,
+        N,
+        datato,
+        Nxto,Nyto,
+        datafrom,
+        Nxfrom,Nyfrom,
+        offsetx,offsety
+    );
+
+
+    return;
+}
+
+void amsfloatimage_region_castcopy1_tf(
+    int threadnum,
+    int nthreads,
+    float *datato,
+    int Nxto,
+    int Nyto,
+    const uint8_t *datafrom,
+    int Nxfrom,
+    int Nyfrom,
+    int offsetx,
+    int offsety
+)
+{
+    int64_t I,I0,I1,Is,N,Ix,Iy,Ia,Ib;
+    int dx,dy;
+
+    float r,g,b,a;
+
+    dx = Nxfrom;
+    dx = (dx>(Nxto-offsetx)) ? (Nxto-offsetx) : dx;
+    dx = (dx<0) ? 0 : dx;
+
+    dy = Nyfrom;
+    dy = (dy>(Nyto-offsety)) ? (Nyto-offsety) : dy;
+    dy = (dy<0) ? 0 : dy;
+
+    N = dx*dy;
+
+    Is = N/nthreads; Is = (Is<1) ? 1 : N;
+    I0 = (threadnum)*Is;
+    I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+    for(I=I0;I<I1;I++)
+    {
+        Ix = I%dx;
+        Iy = I/dx;
+
+        Ia = (Ix+offsetx) + Nxto*(Iy+offsety);
+        Ib =  Ix + Nxfrom*Iy;
+
+        r = (float)(datafrom[0 + 4*Ib])/255.0;
+        g = (float)(datafrom[1 + 4*Ib])/255.0;
+        b = (float)(datafrom[2 + 4*Ib])/255.0;
+        a = (float)(datafrom[3 + 4*Ib])/255.0;
+
+        datato[0 + 4*Ia] = r;
+        datato[1 + 4*Ia] = g;
+        datato[2 + 4*Ia] = b;
+        datato[3 + 4*Ia] = a;
+    }
+
+    return;
+}
+
+void amsfloatimage_region_castcopy2_tf(
+    int threadnum,
+    int nthreads,
+    uint8_t *datato,
+    int Nxto,
+    int Nyto,
+    const float *datafrom,
+    int Nxfrom,
+    int Nyfrom,
+    int offsetx,
+    int offsety
+)
+{
+    int64_t I,I0,I1,Is,N,Ix,Iy,Ia,Ib;
+    int dx,dy;
+
+    int ri,gi,bi,ai;
+
+    dx = Nxfrom;
+    dx = (dx>(Nxto-offsetx)) ? (Nxto-offsetx) : dx;
+    dx = (dx<0) ? 0 : dx;
+
+    dy = Nyfrom;
+    dy = (dy>(Nyto-offsety)) ? (Nyto-offsety) : dy;
+    dy = (dy<0) ? 0 : dy;
+
+    N = dx*dy;
+
+    Is = N/nthreads; Is = (Is<1) ? 1 : N;
+    I0 = (threadnum)*Is;
+    I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+    for(I=I0;I<I1;I++)
+    {
+        Ix = I%dx;
+        Iy = I/dx;
+
+        Ia = (Ix+offsetx) + Nxto*(Iy+offsety);
+        Ib =  Ix + Nxfrom*Iy;
+
+        ri = (int)(datafrom[0 + 4*Ib]*255.0);
+        gi = (int)(datafrom[1 + 4*Ib]*255.0);
+        bi = (int)(datafrom[2 + 4*Ib]*255.0);
+        ai = (int)(datafrom[3 + 4*Ib]*255.0);
+        ri = (ri<0) ? 0 : ri;
+        gi = (ri<0) ? 0 : gi;
+        bi = (ri<0) ? 0 : bi;
+        ai = (ri<0) ? 0 : ai;
+        ri = (ri>255) ? 255 : ri;
+        gi = (ri>255) ? 255 : gi;
+        bi = (ri>255) ? 255 : bi;
+        ai = (ri>255) ? 255 : ai;
+
+        datato[0 + 4*Ia] = (uint8_t)ri;
+        datato[1 + 4*Ia] = (uint8_t)gi;
+        datato[2 + 4*Ia] = (uint8_t)bi;
+        datato[3 + 4*Ia] = (uint8_t)ai;
+    }
+
+    return;
+}
+
+void amsfloatimage_region_castcopy(
+    float *datato,
+    int Nxto,
+    int Nyto,
+    const uint8_t *datafrom,
+    int Nxfrom,
+    int Nyfrom,
+    int offsetx,
+    int offsety
+)
+{
+    int dx,dy;
+    int64_t N;
+
+    if(datato==NULL) return;
+    if(datafrom==NULL) return;
+
+    dx = Nxfrom;
+    dx = (dx>(Nxto-offsetx)) ? (Nxto-offsetx) : dx;
+    dx = (dx<0) ? 0 : dx;
+
+    dy = Nyfrom;
+    dy = (dy>(Nyto-offsety)) ? (Nyto-offsety) : dy;
+    dy = (dy<0) ? 0 : dy;
+
+    N = dx*dy;
+
+    threaded_execute(
+        amsfloatimage_region_castcopy1_tf,
+        N,
+        datato,
+        Nxto,Nyto,
+        datafrom,
+        Nxfrom,Nyfrom,
+        offsetx,offsety
+    );
+
+
+    return;
+}
+
+void amsfloatimage_region_castcopy(
+    uint8_t *datato,
+    int Nxto,
+    int Nyto,
+    const float *datafrom,
+    int Nxfrom,
+    int Nyfrom,
+    int offsetx,
+    int offsety
+)
+{
+    int dx,dy;
+    int64_t N;
+
+    if(datato==NULL) return;
+    if(datafrom==NULL) return;
+
+    dx = Nxfrom;
+    dx = (dx>(Nxto-offsetx)) ? (Nxto-offsetx) : dx;
+    dx = (dx<0) ? 0 : dx;
+
+    dy = Nyfrom;
+    dy = (dy>(Nyto-offsety)) ? (Nyto-offsety) : dy;
+    dy = (dy<0) ? 0 : dy;
+
+    N = dx*dy;
+
+    threaded_execute(
+        amsfloatimage_region_castcopy2_tf,
+        N,
+        datato,
+        Nxto,Nyto,
+        datafrom,
+        Nxfrom,Nyfrom,
+        offsetx,offsety
+    );
+
+
+    return;
+}
+
+void amsfloatimage_region_set_tf(
+    int threadnum,
+    int nthreads,
+    float *data,
+    int Nx, int Ny,
+    int x0, int y0,
+    int x1, int y1,
+    amsfloatpixel val
+)
+{
+    int64_t I,I0,I1,Is,N,Ix,Iy;
+    int dx,dy;
+
+    dx = (x1-x0); dy = (y1-y0);
+    dx = (dx<0) ? 0 : dx;
+    dy = (dy<0) ? 0 : dy;
+    N = dx*dy;
+
+    Is = N/nthreads; Is = (Is<1) ? 1 : N;
+    I0 = (threadnum)*Is;
+    I1 = (threadnum<nthreads-1) ? (threadnum+1)*Is : N;
+
+    for(I=I0;I<I1;I++)
+    {
+        Ix = I%dx;
+        Iy = I/dx;
+
+        data[0 + 4*(Ix+x0) + 4*Nx*(Iy+y0)] = val.R;
+        data[1 + 4*(Ix+x0) + 4*Nx*(Iy+y0)] = val.G;
+        data[2 + 4*(Ix+x0) + 4*Nx*(Iy+y0)] = val.B;
+        data[3 + 4*(Ix+x0) + 4*Nx*(Iy+y0)] = val.A;
+    }
+
+    return;
+}
+
+void amsfloatimage_region_set(
+    float *data,
+    int Nx, int Ny,
+    int x0, int y0,
+    int x1, int y1,
+    amsfloatpixel val
+)
+{
+    int dx,dy;
+    int64_t N;
+
+    if(data==NULL) return;
+
+    dx = (x1-x0); dy = (y1-y0);
+    dx = (dx<0) ? 0 : dx;
+    dy = (dy<0) ? 0 : dy;
+    N = dx*dy;
+
+    threaded_execute(
+        amsfloatimage_region_set_tf,
+        N,
+        data, Nx,Ny, 
+        x0, y0, x1, y1, 
+        val
+    );
+
+    return;
+}
+
+
 }; //end namespace imglib4
 }; //end namespace ams
\ No newline at end of file