dolfin team mailing list archive

[A Navaei <axnavaei@xxxxxxxxxxxxxx>]

I applied a slight correction to the inline c++ code to get the input
image mapped correctly to the function (the indices were not in the
right order) -- see the attachment.

For a complete c++ implementation, instead of sub-classing Function, I
tried sub-classing DiscreteFunction and using it in Function so that
the image-to-function class does not return a user defined class.
However, I was not able to find a pure c++ example (without using
Forms) demonstrating how to create a Function from scratch, specially,
how to create a Dofmap, or avoid one if we do not need it. Here is
where a Dofmap instance is required:

http://code.google.com/p/wrapitk/source/browse/trunk/ExternalProjects/ItkDolfin/src/itkImageToDolfinFunction.txx?r=303#13

ImageDiscreteFunction is defined here:

http://code.google.com/p/wrapitk/source/browse/trunk/ExternalProjects/ItkDolfin/src/itkImageToDolfinFunction.h?r=303#83

This interface is able to convert itk::Image to dolfin::Function and
work for 2D and 3D data. from the itk side, any nD images can be
introduced, but I don't think if we are interested in any other
dimensions in dolfin? Also, it will support the conversion of vector
images. The reverse conversion will be added soon.

-Ali

2009/2/7  <kent-and@xxxxxxxxx>:
> Here is the same example using a compiled function.
>
> Kent
>
>> Do you have python-itk-numpy installed?
>>
>> 2009/2/6 Kent Andre <kent-and@xxxxxxxxx>:
>>>
>>> Great, but I get the following problem:
>>>
>>>
>>>
>>> Traceback (most recent call last):
>>>  File "itk-dolfin.py", line 16, in <module>
>>>    itk2numpy = itk.PyBuffer[inType]
>>>  File "/usr/lib/InsightToolkit/WrapITK/Python/itkLazy.py", line 14, in
>>> __getattribute__
>>>    value = types.ModuleType.__getattribute__(self, attr)
>>> AttributeError: 'LazyITKModule' object has no attribute 'PyBuffer'
>>>
>>> On to., 2009-02-05 at 19:18 +0000, A Navaei wrote:
>>>> >
>>>> > Let me know when you have created debian packages. It seems I have to
>>>> > rebuild everything, even cmake. Compilation is in principle no
>>>> problem,
>>>> > but I already compile a decent number of packages on a regular basis
>>>> and
>>>> > being non-bleeding edge on packages I don't know well seems tempting.
>>>> >
>>>> > I might download and compile, but not today.
>>>> >
>>>> > Kent
>>>> >
>>>> >
>>>>
>>>> python-itk debian, based on the new WrapITK, is just published:
>>>> http://code.google.com/p/wrapitk/wiki/WrapITKBinaries .The 32 bit
>>>> version is ready, the 64 bit version is being built on opensuse build
>>>> service and will be available in a few hours.
>>>>
>>>> I have tried the [itk]->[numpy]->[dolfin] example, and it already
>>>> looks great, here is the output:
>>>>
>>>> Checking mesh ordering (finished).
>>>> 211 * 212 = 44732 pixels took 246.780 ms transfering from itk to
>>>> dolfin through numpy
>>>> That is 5.516854 us per pixel.
>>>>
>>>> I am going to add the itk-dolfin interface in c++ in the next few
>>>> days, which should make the conversion even faster leading to more
>>>> reasonable performance for 3D data. Below is an updated version of the
>>>> example. Let me know if there are any problems.
>>>>
>>>>
>>>> -Ali
>>>>
>>>> ----------------------------------------------------------
>>>> # [itk]->[numpy]->[dolfin] test
>>>>
>>>> from dolfin import *
>>>> from numpy import *
>>>> import itk
>>>> import time
>>>>
>>>> dim = 2
>>>> inType = itk.Image[itk.D, dim]
>>>>
>>>> reader = itk.ImageFileReader[inType].New(FileName="/path/to/file.bmp")
>>>> reader.Update()
>>>>
>>>> t1 = time.time()
>>>>
>>>> itk2numpy = itk.PyBuffer[inType]
>>>> numpy_arr = itk2numpy.GetArrayFromImage( reader.GetOutput() )
>>>> shape = numpy_arr.shape
>>>>
>>>> mesh = UnitSquare(shape[0]-1, shape[1]-1)
>>>>
>>>> class ImageFunction(Function):
>>>>     def eval(self, value, x):
>>>>         i = int((self.A.shape[0]-1)*x[0])
>>>>         j = int((self.A.shape[1]-1)*x[1])
>>>> #        print i,j
>>>>         value[0] = self.A[(i,j)]
>>>>
>>>> V = FunctionSpace(mesh, "CG", 1)
>>>> f = ImageFunction(V)
>>>> f.A = numpy_arr
>>>>
>>>> t2 = time.time()
>>>> print '%i * %i = %i pixels took %0.3f ms transfering from itk to
>>>> dolfin through numpy'% (shape[0], shape[1], shape[0]*shape[1],
>>>> (t2-t1)*1e3)
>>>> print 'That is %f us per pixel.'% ((t2-t1) * 1e6 / (shape[0]*shape[1]))
>>>>
>>>> plot(f) # can viper plot 2d images better than this?
>>>> interactive()
>>>
>>>
>>
>

from dolfin import *
from numpy import *
import itk
import time

dim = 2
inType = itk.Image[itk.D, dim]

reader = itk.ImageFileReader[inType].New(FileName="/path/to/image.bmp")
reader.Update()

t1 = time.time()

itk2numpy = itk.PyBuffer[inType]
numpy_arr = itk2numpy.GetArrayFromImage( reader.GetOutput() )
shape = numpy_arr.shape

mesh = UnitSquare(shape[0]-1, shape[1]-1)

code = """
class ItkFunc : public Function
{
    private:
        int N; 
        int M;
        double *image; 

    public: 

        void set_image(PyObject* a_){ 
            // Do some more type checks ... 
            PyArrayObject* a=(PyArrayObject*) a_;
            N = a->dimensions[0];
            M = a->dimensions[1];
            image = (double*) a->data; 
        }



        ItkFunc(FunctionSpace& V) : Function(V), N(0), M(0) {}

        void eval(double* values, const double* x) const
        {
            int i = int((M-1)*x[0]);
            int j = int((N-1)*x[1]);
            int ij = i + (M)*j;
            values[0] = image[ij]; 
        }

};
"""

V = FunctionSpace(mesh, "CG", 1)
f = Function(V, code)
f.set_image(numpy_arr)

t2 = time.time()
print '%i * %i = %i pixels took %0.3f ms transfering from itk to dolfin through numpy'% (shape[0], shape[1], shape[0]*shape[1], (t2-t1)*1e3)
print 'That is %f us per pixel.'% ((t2-t1) * 1e6 / (shape[0]*shape[1]))

plot(f) # can viper plot 2d images better than this?
interactive()