Detailed Description

The image registration problem is modeled here with the finite element method. Image registration is, in general, an ill-posed problem. Thus, we use an optimization scheme where the optimization criterion is given by a regularized variational energy. The variational energy arises from modeling the image as a physical body on which external forces act. The body is allowed to deform so as to minimize the applied force. The resistance of the physical body to deformation, determined by the physics associated with the body, serves to regularize the solution. The forces applied to the body are, generally, highly non-linear and so the body is allowed to deform slowly and incrementally. The direction it deforms follows the gradient of the potential energy (the force) we define. The potential energies we may choose from are given by the itk image-to-image metrics. The choices and the associated direction of descent are : Mean Squares (minimize), Normalized Cross-Correlation (maximize) Mutual Information (maximize). Note that we have to set the direction (SetDescentDirection) when we choose a metric. The forces driving the problem may also be given by user-supplied landmarks. The corners of the image, in this example, are always pinned. This example is designed for 2D or 3D images. A rectilinear mesh is generated automatically given the correct element type (Quadrilateral or Hexahedral).

Our specific Solver for this example uses trapezoidal time stepping. This is a method for solving a second-order PDE in time. The solution is penalized by the zeroth (mass matrix) and first derivatives (stiffness matrix) of the shape functions. There is an option to perform a line search on the energy after each iteration. Optimal parameter settings require experimentation. The following approach tends to work well : Choose the relative size of density to elasticity (e.g. Rho / E ~= 1.) such that the image deforms locally and slowly. This also affects the stability of the solution. Choose the time step to control the size of the deformation at each step. Choose enough iterations to allow the solution to converge (this may be automated).

Reading images is up to the user. Either set the images using SetMoving/FixedImage or see the ReadImages function.

Note:

This code works for only 2 or 3 dimensions b/c we do not have > 3D elements.

TODO : Keep the full field around (if using re-gridding). Introduce compensation for kinematic non-linearity in time (if using Eulerian frame).

Allow people to add/remove/invoke observers (callbacks) to any ITK object. This is an implementation of the subject/observer design pattern. An observer is added by specifying an event to respond to and an itk::Command to execute. It returns an unsigned long tag which can be used later to remove the event or retrieve the command. The memory for the Command becomes the responsibility of this object, so don't pass the same instance of a command to two different objects

The GenerateData method normally allocates the buffers for all of the outputs of a filter. Some filters may want to override this default behavior. For example, a filter may have multiple outputs with varying resolution. Or a filter may want to process data in place by grafting its input to its output.

This method is called when itkExceptionMacro executes. It allows the debugger to break on error.

virtual void itk::ProcessObject::CacheInputReleaseDataFlags

(

)

[protected, virtual, inherited]

Cache the state of any ReleaseDataFlag's on the inputs. While the filter is executing, we need to set the ReleaseDataFlag's on the inputs to false in case the current filter is implemented using a mini-pipeline (which will try to release the inputs). After the filter finishes, we restore the state of the ReleaseDataFlag's before the call to ReleaseInputs().

This function calls the actual region copier to do the mapping from input image space to output image space. It uses a Function object used for dispatching to various routines to copy an input region (start index and size) to an output region. For most filters, this is a trivial copy because most filters require the input dimension to match the output dimension. However, some filters like itk::UnaryFunctorImageFilter can support output images of a higher dimension that the input.

The default copier uses a "dispatch pattern" to call one of three overloaded functions depending on whether the input and output images are the same dimension, the input is a higher dimension that the output, or the input is of a lower dimension than the output. The use of an overloaded function is required for proper compilation of the various cases.

For the latter two cases, trivial implementations are used. If the input image is a higher dimension than the output, the first portion of the input region is copied to the output region. If the input region is a lower dimension than the output, the input region information is copied into the first portion of the output region and the rest of the output region is set to zero.

If a filter needs a different default behavior, it can override this method.

This function calls the actual region copier to do the mapping from output image space to input image space. It uses a Function object used for dispatching to various routines to copy an output region (start index and size) to an input region. For most filters, this is a trivial copy because most filters require the input dimension to match the output dimension. However, some filters like itk::ExtractImageFilter can support output images of a lower dimension that the input.

This function object can be used by GenerateOutputInformation() to copy the input LargestPossibleRegion to the output LargestPossibleRegion and can also be used in GenerateData or ThreadedGenerateData() where a filter may need to map an output region to an input region.

The default copier uses a "dispatch pattern" to call one of three overloaded functions depending on whether the input and output images are the same dimension, the input is a higher dimension that the output, or the input is of a lower dimension than the output. The use of an overloaded function is required for proper compilation of the various cases.

For the latter two cases, trivial implementations are used. If the input image is a higher dimension than the output, the output region information is copied into the first portion of the input region and the rest of the input region is set to zero. If the input region is a lower dimension than the output, the first portion of the output region is copied to the input region.

If a filter needs a different default behavior, it can override this method. The ExtractImageFilter overrides this function object so that if the input image is a higher dimension than the output image, the filter can control "where" in the input image the output subimage is extracted (as opposed to mapping to first few dimensions of the input).

Create an object from an instance, potentially deferring to a factory. This method allows you to create an instance of an object that is exactly the same type as the referring object. This is useful in cases where an object has been cast back to a base class.

We check the jacobian of the current deformation field. If it is < threshold, we begin diffeomorphism enforcement: 1) Warp the moving image. 2) Set the vector field to zero. 3) Set the warped moving image as the new moving image, resizing if necessary.

Give the process object a chance to indictate that it will produce more output than it was requested to produce. For example, many imaging filters must compute the entire output at once or can only produce output in complete slices. Such filters cannot handle smaller requested regions. These filters must provide an implementation of this method, setting the output requested region to the size they will produce. By default, a process object does not modify the size of the output requested region.

A version of GenerateData() specific for image processing filters. This implementation will split the processing across multiple threads. The buffer is allocated by this method. Then the BeforeThreadedGenerateData() method is called (if provided). Then, a series of threads are spawned each calling ThreadedGenerateData(). After all the threads have completed processing, the AfterThreadedGenerateData() method is called (if provided). If an image processing filter cannot be threaded, the filter should provide an implementation of GenerateData(). That implementation is responsible for allocating the output buffer. If a filter an be threaded, it should NOT provide a GenerateData() method but should provide a ThreadedGenerateData() instead.

What is the input requested region that is required to produce the output requested region? The base assumption for image processing filters is that the input requested region can be set to match the output requested region. If a filter requires more input (for instance a filter that uses neighborhoods needs more input than output to avoid introducing artificial boundary conditions) or less input (for instance a magnify filter) will have to override this method. In doing so, it should call its superclass' implementation as its first step. Note that imaging filters operate differently than the classes to this point in the class hierachy. Up till now, the base assumption has been that the largest possible region will be requested of the input.

Generate the information decribing the output data. The default implementation of this method will copy information from the input to the output. A filter may override this method if its output will have different information than its input. For instance, a filter that shrinks an image will need to provide an implementation for this method that changes the spacing of the pixels. Such filters should call their superclass' implementation of this method prior to changing the information values they need (i.e. GenerateOutputInformation() should call Superclass::GenerateOutputInformation() prior to changing the information.

Given one output whose requested region has been set, how should the requested regions for the remaining outputs of the process object be set? By default, all the outputs are set to the same requested region. If a filter needs to produce different requested regions for each output, for instance an image processing filter producing several outputs at different resolutions, then that filter may override this method and set the requested regions appropriatedly.

Note that a filter producing multiple outputs of different types is required to override this method. The default implementation can only correctly handle multiple outputs of the same type.

Get the command associated with the given tag. NOTE: This returns a pointer to a Command, but it is safe to asign this to a Command::Pointer. Since Command inherits from LightObject, at this point in the code, only a pointer or a reference to the Command can be used.

Get the size of the input vector. This is merely the size of the input vector, not the number of inputs that have valid DataObject's assigned. Use GetNumberOfValidRequiredInputs() to determine how many inputs are non-null.

Get the number of valid inputs. This is the number of non-null entries in the input vector in the first NumberOfRequiredInputs slots. This method is used to determine whether the necessary required inputs have been set. Subclasses of ProcessObject may override this implementation if the required inputs are not the first slots in input vector.

Get the output data of this process object. The output of this function is not valid until an appropriate Update() method has been called, either explicitly or implicitly. Both the filter itself and the data object have Update() methods, and both methods update the data. Here are three ways to use GetOutput() and make sure the data is valid. In these examples, image is a pointer to some Image object, and the particular ProcessObjects involved are filters. The same examples apply to non-image (e.g. Mesh) data as well.

In this situation, someFilter and anotherFilter are said to constitute a pipeline.

image = someFilter->GetOutput();
image->Update();

someFilter->Update();
image = someFilter->GetOutput();

(In the above example, the two lines of code can be in either order.)

Note that Update() is not called automatically except within a pipeline as in the first example. When streaming (using a StreamingImageFilter) is activated, it may be more efficient to use a pipeline than to call Update() once for each filter in turn.

Graft the specified data object onto this ProcessObject's idx'th output. This is similar to the GraftOutput method except it allows you to specify which output is affected. The specified index must be a valid output number (less than ProcessObject::GetNumberOfOutputs()). See the GraftOutput for general usage information.

Graft the specified DataObject onto this ProcessObject's output. This method grabs a handle to the specified DataObject's bulk data to used as its output's own bulk data. It also copies the region ivars (RequestedRegion, BufferedRegion, LargestPossibleRegion) and meta-data (Spacing, Origin) from the specified data object into this filter's output data object. Most importantly, however, it leaves the Source ivar untouched so the original pipeline routing is intact. This method is used when a process object is implemented using a mini-pipeline which is defined in its GenerateData() method. The usage is:

// setup the mini-pipeline to process the input to this filter
firstFilterInMiniPipeline->SetInput( this->GetInput() );
// setup the mini-pipeline to calculate the correct regions// and write to the appropriate bulk data block
lastFilterInMiniPipeline->GraftOutput( this->GetOutput() );
// execute the mini-pipeline
lastFilterInMiniPipeline->Update();
// graft the mini-pipeline output back onto this filter's output.// this is needed to get the appropriate regions passed back.
this->GraftOutput( lastFilterInMiniPipeline->GetOutput() );

Methods invoked by Print() to print information about the object including superclasses. Typically not called by the user (use Print() instead) but used in the hierarchical print process to combine the output of several classes.

PushBackInput(), PushFronInput() in the public section force the input to be the type expected by an ImageToImageFilter. However, these methods end of "hiding" the versions from the superclass (ProcessObject) whose arguments are DataObjects. Here, we re-expose the versions from ProcessObject to avoid warnings about hiding methods from the superclass.

Push/Pop the input of this process object. These methods allow a filter to model its input vector as a queue or stack. These routines may not be appropriate for all filters, especially filters with different types of inputs. These routines follow the semantics of STL.

The routines are useful for applications that need to process "rolling" sets of images. For instance, if an application has 10 images and they need to run a filter on images 1, 2, 3, 4, then run the filter on images 2, 3, 4, 5, then run the filter on images 3, 4, 5, 6, the application can accomplish this by popping an input off the front of the input list and push a new image onto the back of input list. Again, this only makes sense for filters that single type of input.

Other uses are also possible. For a single input filter, pushing and popping inputs allow the application to temporarily replace an input to a filter.

PushBackInput(), PushFronInput() in the public section force the input to be the type expected by an ImageToImageFilter. However, these methods end of "hiding" the versions from the superclass (ProcessObject) whose arguments are DataObjects. Here, we re-expose the versions from ProcessObject to avoid warnings about hiding methods from the superclass.

A filter may need to release its input's bulk data after it has finished calculating a new output. The filter may need to release the inputs because the user has turned on the ReleaseDataFlag or it may need to release the inputs because the filter is an "in place" filter and it has overwritten its input with its output data. The implementation here simply checks the ReleaseDataFlag of the inputs. InPlaceImageFilter overrides this method so release the input it has overwritten.

The FEM filter can generate its own mesh for 2 or 3 dimensions, if none is provided. The mesh is generated for quadrilaterals in 2D and hexahedra in 3D. This function sets the number of elements generated along each dimension at the resolution designated by "which". E.g. to generate 10 pixels per element in each dimension in the 1st resolution, use SetMeshResolution(10,0);.

This determines the number of integration points to use at each resolution. These integration points are used to generate the force. The actual number used will be i^d, where d is the number of parameters in the elements local domain.

Set the execution progress of a process object. The progress is a floating number in [0,1] with 0 meaning no progress and 1 meaning the filter has completed execution. The ProgressEvent is NOT invoked.

Turn on/off the flags to control whether the bulk data belonging to the outputs of this ProcessObject are released/reallocated during an Update(). In limited memory scenarios, a user may want to force the elements of a pipeline to release any bulk data that is going to be regenerated anyway during an Update() in order to control peak memory allocation. Note that this flag is different from the ReleaseDataFlag. ReleaseDataFlag manages the deallocation of a ProcessObject's bulk output data once that data has been consumed by a downstream ProcessObject. The ReleaseDataBeforeUpdateFlag manages the deallocation/reallocation of bulk data during a pipeline update to control peak memory utilization. Default value is on.

virtual void itk::ProcessObject::SetReleaseDataFlag

(

bool

flag

)

[virtual, inherited]

Turn on/off the flags to control whether the bulk data belonging to the outputs of this ProcessObject are released after being used by a downstream ProcessObject. Default value is off. Another options for controlling memory utilization is the ReleaseDataBeforeUpdateFlag.

The warped reference image will be written to this file name with the extension "11.img" appended to it. One can also output the image after every iteration, yielding result11.img, result12.img, etc. by uncommenting the code at the end of IterativeSolve.

Split the output's RequestedRegion into "num" pieces, returning region "i" as "splitRegion". This method is called "num" times. The regions must not overlap. The method returns the number of pieces that the routine is capable of splitting the output RequestedRegion, i.e. return value is less than or equal to "num".

If an imaging filter can be implemented as a multithreaded algorithm, the filter will provide an implementation of ThreadedGenerateData(). This superclass will automatically split the output image into a number of pieces, spawn multiple threads, and call ThreadedGenerateData() in each thread. Prior to spawning threads, the BeforeThreadedGenerateData() method is called. After all the threads have completed, the AfterThreadedGenerateData() method is called. If an image processing filter cannot support threading, that filter should provide an implementation of the GenerateData() method instead of providing an implementation of ThreadedGenerateData(). If a filter provides a GenerateData() method as its implementation, then the filter is responsible for allocating the output data. If a filter provides a ThreadedGenerateData() method as its implementation, then the output memory will allocated automatically by this superclass. The ThreadedGenerateData() method should only produce the output specified by "outputThreadRegion" parameter. ThreadedGenerateData() cannot write to any other portion of the output image (as this is responsibility of a different thread).

Bring this filter up-to-date. Update() checks modified times against last execution times, and re-executes objects if necessary. A side effect of this method is that the whole pipeline may execute in order to bring this filter up-to-date. This method updates the currently prescribed requested region. If no requested region has been set on the output, then the requested region will be set to the largest possible region. Once the requested region is set, Update() will make sure the specified requested region is up-to-date. This is a confusing side effect to users who are just calling Update() on a filter. A first call to Update() will cause the largest possible region to be updated. A second call to Update() will update that same region. If a modification to the upstream pipeline cause a filter to have a different largest possible region, this second call to Update() will not cause the output requested region to be reset to the new largest possible region. Instead, the output requested region will be the same as the last time Update() was called. To have a filter always to produce its largest possible region, users should call UpdateLargestPossibleRegion() instead.

Like Update(), but sets the output requested region to the largest possible region for the output. This is the method users should call if they want the entire dataset to be processed. If a user wants to update the same output region as a previous call to Update() or a previous call to UpdateLargestPossibleRegion(), then they should call the method Update().

Update the information decribing the output data. This method transverses up the pipeline gathering modified time information. On the way back down the pipeline, this method calls GenerateOutputInformation() to set any necessary information about the output data objects. For instance, a filter that shrinks an image will need to provide an implementation for GenerateOutputInformation() that changes the spacing of the pixels. Such filters should call their superclass' implementation of GenerateOutputInformation prior to changing the information values they need (i.e. GenerateOutputInformation() should call Superclass::GenerateOutputInformation() prior to changing the information.