dolfin team mailing list archive

[kent-and@xxxxxxxxx]

> On Mon, Jan 26, 2009 at 03:33:09PM +0000, A Navaei wrote:
>> Anders,
>>
>> Kent and I have discussed connecting itk and dolfin and, following
>> your advice, one viable way of itk->dolfin data flow can be summarised
>> in the following c++ snippet:
>>
>> -------------------------------------
>> typedef double PixelType
>> const unsigned int Dim = 3;
>> typedef itk::Image< PixelType, Dim > ImageType;
>>
>> class F : public Function
>> {
>> public:
>>  SetImageData(ImageType id) {this->ImageData = id;};
>>
>>  void eval(double* values, const double* x) const
>>  {
>>    values[0] = ImageData.GetPixel(x);
>>  }
>>
>> protected:
>>  ImageType ImageData;
>> };
>> --------------------------------------
>>
>> The problem with the above approach is that direct access of pixel
>> data with GetPixel() in itk, and generally in a generic image data
>> structure, is not efficient. Typically, in a 2D or 3D image, this
>> method is, correspondingly, called 10e4 to 10e6 times leading to slow
>> data structure conversion. This is even more pronounces when the
>> conversion has to happen many times.
>>
>> There are at least 2 alternatives for this:
>>
>> (1) Passing the image data pointer to a Function class.
>>
>> (2) Looping in the itk API rather than the dolfin API. The c++
>> iterator implemented in itk provides a much more efficient interface
>> for data assignment, however, we need to have access to dolfin
>> Function data.
>>
>> In both of the above approaches, we assume that the image data has a
>> complete rectangular shape with known sizes which provides a simple
>> look-up rule. Is it currently possible to implement either of the
>> above methods? If not, how can dolfin be extended to support this?
>
> I suggest just copying all image data to an array in the constructor
> so you just access it once.
>
> Then in your eval function you just lookup the values in the array.
> Depending on how you represent your function, you might want to
> overload the following eval function:
>
>   void eval(double* values, const Data& data) const
>   {
>
>   }
>
> You can access the current cell by data.cell(). If you have a DG(0)
> representation of your image field, you may then simply map the cell
> indices to the pixel values.
>
> --
> Anders

I agree with Anders.

The dof_map does of course have a very simple structure in this case, which
can be written as an explicit formula. This is an advantage we usually don't
have when working on unstructured meshes. Still, I think the cpu time and
memory
will not be dominated by the dof_map so it might be smart to wait about
taking
advantage of this until after profiling.

Can you store the images eg as a numpy array?

Kent