[vtkusers] double vtkDataArray::GetTuple1(vtkIdType i) vs void vtkDataArray::GetTuple(vtkIdType i, double *tuple)

[David Gobbi]

Hi Maarten,

I actually did write an iterator with such a function, but I never got
around to contributing it to VTK.  I called it vtkImageRegionIterator.
The source code is here:
http://github.com/dgobbi/AIRS/tree/master/ImageSegmentation

It has three separate inline methods for getting the voxel:
GetIndexX(), GetIndexY(), GetIndexZ().

A couple years ago I had some big plans for improving the VTK
image iterators, but like so many things, there was never enough
time.

  David



On Thu, May 8, 2014 at 7:20 AM, Maarten Beek <beekmaarten at yahoo.com> wrote:
> Hi David,
>
> Thanks so much for your clarification!
> I did see issues in the resulting image but was attributing them to noise in
> the original images. Now I see my code just had a huge bug in it.
>
> I guess I would need a row counter in order to initialize voxelid correctly.
>
> Wouldn't it be nice if the image iterator had a function like getVoxel(int
> ijk[3])...? I believe ITK has such an iterator.
>
> Maarten
>
>
> On Wednesday, May 7, 2014 5:26:23 PM, David Gobbi <david.gobbi at gmail.com>
> wrote:
> I should add: you actually have to re-initialize the voxelid at the
> beginning of each row (i.e. at the beginning of each span of the
> iterator).
>
> On Wed, May 7, 2014 at 3:22 PM, David Gobbi <david.gobbi at gmail.com> wrote:
>> Hi Maarten,
>>
>> Since you know it's a vtkUnsignedCharArray, you shouldn't be using
>> the SetTuple/GetTuple methods.  Use the SetTupleValue and
>> GetTupleValue methods, they're much faster because 1) they aren't
>> virtual methods and 2) they don't have to convert to/from double.
>> Or even better, if you know that your array only has one component,
>> you can use SetValue/GetValue.
>>
>> Your math for voxelid is too simplistic.  Depending on how the update
>> extent is divided amongst the threads, it might give you a location
>> that is totally wrong.  To see why, draw a square on a piece of paper.
>> That square represents the total extent of the output.  Then divide
>> that square into quadrants, and imagine that each quadrant goes
>> to a different Execute method.  Then think about how each thread
>> is going to march through the data.
>>
>>  David
>>
>>
>> On Wed, May 7, 2014 at 2:20 PM, Maarten Beek <beekmaarten at yahoo.com>
>> wrote:
>>> Hi all,
>>>
>>> I experienced a crash in my algorithm derived from
>>> vtkThreadedImageAlgorithm, that (I think) I solved, but that I don't
>>> understand. Is there some one who can shed some light on the following?
>>>
>>> What I am trying to do is adding a flag identifying the input(s) that
>>> contribute to the value of each pixel in the output image
>>>
>>>
>>> //----------------------------------------------------------------------------
>>> // This templated function executes the filter for any type of data.
>>> template <class OT>
>>> void cmvtkTemplatedExecute(cmvtkDRRGenerator2* self,
>>>                                vtkDataSetCollection* inputs,
>>> vtkImageData*
>>> output,
>>>                                int extent[6], int threadId, OT*)
>>> {
>>>    vtkUnsignedCharArray* uchararray =
>>>
>>> vtkUnsignedCharArray::SafeDownCast(output->GetPointData()->GetScalars("input_flag"));
>>>
>>>    int numinputs = inputs->GetNumberOfItems();
>>>    for( int input_i = 0; input_i < numinputs; ++input_i )
>>>    {
>>>        vtkImageData* input =
>>> vtkImageData::SafeDownCast(inputs->GetItem(input_i));
>>>        if( 0 == input )
>>>        {
>>>            // should not happen
>>>            continue;
>>>        }
>>>        const int* dims= output->GetDimensions();
>>>
>>>        vtkImageProgressIterator<OT> outIt(output, extent, self,
>>> threadId);
>>>        vtkIdType voxelid = extent[0] + size[0]*extent[2] +
>>> size[0]*size[1]*extent[4];
>>>
>>>        // Loop through output pixels
>>>        while( !outIt.IsAtEnd() )
>>>        {
>>>            OT* outSI = outIt.BeginSpan();
>>>            OT* outSIEnd = outIt.EndSpan();
>>>
>>>            while( outSI != outSIEnd )
>>>            {
>>>                // start with setting pixel to 0.0
>>>                if( input_i == 0 )
>>>                {
>>>                    *outSI = 0.0;
>>>                }
>>>
>>>                OT result = 0.0;
>>>
>>>                //.... some math resulting in a value for result....
>>>
>>>                *outSI += result;
>>>
>>>                // this input image participates in voxel values; adjust
>>> bit
>>>                if( 0.0 < result && 0 != uchararray )
>>>                {
>>>                    // The commented code crashes...
>>>
>>>                    //unsigned char oldvalue =
>>> uchararray->GetTuple1(voxelid);
>>>                    //unsigned char newvalue = oldvalue | (1<<input_i);
>>>                    //uchararray->SetTuple1(voxelid, newvalue);
>>>
>>>                    double dOldValue;
>>>                    uchararray->GetTuple(voxelid, &dOldValue);
>>>                    unsigned char chOldValue = static_cast<unsigned
>>> char>(dOldValue);
>>>                    unsigned char chNewValue = chOldValue | (1<<input_i);
>>>                    double dNewValue = static_cast<double>(chNewValue);
>>>                    uchararray->SetTuple(voxelid, &dNewValue);
>>>                }
>>>
>>>                ++outSI;
>>>                ++voxelid;
>>>            }
>>>            outIt.NextSpan();
>>>        }
>>>    } // input_i < numinputs
>>> }
>
>