diff --git a/Core/Code/DataManagement/mitkPlaneGeometry.cpp b/Core/Code/DataManagement/mitkPlaneGeometry.cpp index 290e79a61b..b254aca42e 100644 --- a/Core/Code/DataManagement/mitkPlaneGeometry.cpp +++ b/Core/Code/DataManagement/mitkPlaneGeometry.cpp @@ -1,753 +1,755 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision$ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkPlaneGeometry.h" #include "mitkPlaneOperation.h" #include "mitkInteractionConst.h" #include "mitkLine.h" #include #include namespace mitk { mitk::PlaneGeometry::PlaneGeometry() { Initialize(); } mitk::PlaneGeometry::~PlaneGeometry() { } void PlaneGeometry::Initialize() { Superclass::Initialize(); } void PlaneGeometry::EnsurePerpendicularNormal(mitk::AffineTransform3D *transform) { //ensure row(2) of transform to be perpendicular to plane, keep length. VnlVector normal = vnl_cross_3d( transform->GetMatrix().GetVnlMatrix().get_column(0), transform->GetMatrix().GetVnlMatrix().get_column(1) ); normal.normalize(); ScalarType len = transform->GetMatrix() .GetVnlMatrix().get_column(2).two_norm(); if (len==0) len = 1; normal*=len; Matrix3D matrix = transform->GetMatrix(); matrix.GetVnlMatrix().set_column(2, normal); transform->SetMatrix(matrix); } void PlaneGeometry::SetIndexToWorldTransform(mitk::AffineTransform3D *transform) { EnsurePerpendicularNormal(transform); Superclass::SetIndexToWorldTransform(transform); } void PlaneGeometry::SetBounds(const BoundingBox::BoundsArrayType &bounds) { //currently the unit rectangle must be starting at the origin [0,0] assert(bounds[0]==0); assert(bounds[2]==0); //the unit rectangle must be two-dimensional assert(bounds[1]>0); assert(bounds[3]>0); Superclass::SetBounds(bounds); } void PlaneGeometry::IndexToWorld( const Point2D &pt_units, Point2D &pt_mm ) const { pt_mm[0]=m_ScaleFactorMMPerUnitX*pt_units[0]; pt_mm[1]=m_ScaleFactorMMPerUnitY*pt_units[1]; } void PlaneGeometry::WorldToIndex( const Point2D &pt_mm, Point2D &pt_units ) const { pt_units[0]=pt_mm[0]*(1.0/m_ScaleFactorMMPerUnitX); pt_units[1]=pt_mm[1]*(1.0/m_ScaleFactorMMPerUnitY); } void PlaneGeometry::IndexToWorld( const Point2D &atPt2d_units, const Vector2D &vec_units, Vector2D &vec_mm) const { MITK_WARN<<"Warning! Call of the deprecated function PlaneGeometry::IndexToWorld(point, vec, vec). Use PlaneGeometry::IndexToWorld(vec, vec) instead!"; this->IndexToWorld(vec_units, vec_mm); } void PlaneGeometry::IndexToWorld(const Vector2D &vec_units, Vector2D &vec_mm) const { vec_mm[0] = m_ScaleFactorMMPerUnitX * vec_units[0]; vec_mm[1] = m_ScaleFactorMMPerUnitY * vec_units[1]; } void PlaneGeometry::WorldToIndex( const Point2D &atPt2d_mm, const Vector2D &vec_mm, Vector2D &vec_units) const { MITK_WARN<<"Warning! Call of the deprecated function PlaneGeometry::WorldToIndex(point, vec, vec). Use PlaneGeometry::WorldToIndex(vec, vec) instead!"; this->WorldToIndex(vec_mm, vec_units); } void PlaneGeometry::WorldToIndex( const Vector2D &vec_mm, Vector2D &vec_units) const { vec_units[0] = vec_mm[0] * ( 1.0 / m_ScaleFactorMMPerUnitX ); vec_units[1] = vec_mm[1] * ( 1.0 / m_ScaleFactorMMPerUnitY ); } void PlaneGeometry::InitializeStandardPlane( mitk::ScalarType width, ScalarType height, const Vector3D & spacing, PlaneGeometry::PlaneOrientation planeorientation, ScalarType zPosition, bool frontside, bool rotated ) { AffineTransform3D::Pointer transform; transform = AffineTransform3D::New(); AffineTransform3D::MatrixType matrix; AffineTransform3D::MatrixType::InternalMatrixType &vnlmatrix = matrix.GetVnlMatrix(); vnlmatrix.set_identity(); vnlmatrix(0,0) = spacing[0]; vnlmatrix(1,1) = spacing[1]; vnlmatrix(2,2) = spacing[2]; transform->SetIdentity(); transform->SetMatrix(matrix); InitializeStandardPlane(width, height, transform.GetPointer(), planeorientation, zPosition, frontside, rotated); } void PlaneGeometry::InitializeStandardPlane( mitk::ScalarType width, ScalarType height, const AffineTransform3D* transform, PlaneGeometry::PlaneOrientation planeorientation, ScalarType zPosition, bool frontside, bool rotated ) { Superclass::Initialize(); //construct standard view Point3D origin; VnlVector rightDV(3), bottomDV(3); origin.Fill(0); int normalDirection; switch(planeorientation) { case Transversal: if(frontside) { if(rotated==false) { FillVector3D(origin, 0, 0, zPosition); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, 1, 0); } else { FillVector3D(origin, width, height, zPosition); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, -1, 0); } } else { if(rotated==false) { FillVector3D(origin, width, 0, zPosition); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, 1, 0); } else { FillVector3D(origin, 0, height, zPosition); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, -1, 0); } } normalDirection = 2; break; case Frontal: if(frontside) { if(rotated==false) { FillVector3D(origin, 0, zPosition, 0); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, 0, 1); } else { FillVector3D(origin, width, zPosition, height); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, 0, -1); } } else { if(rotated==false) { FillVector3D(origin, width, zPosition, 0); FillVector3D(rightDV, -1, 0, 0); FillVector3D(bottomDV, 0, 0, 1); } else { FillVector3D(origin, 0, zPosition, height); FillVector3D(rightDV, 1, 0, 0); FillVector3D(bottomDV, 0, 0, -1); } } normalDirection = 1; break; case Sagittal: if(frontside) { if(rotated==false) { FillVector3D(origin, zPosition, 0, 0); FillVector3D(rightDV, 0, 1, 0); FillVector3D(bottomDV, 0, 0, 1); } else { FillVector3D(origin, zPosition, width, height); FillVector3D(rightDV, 0, -1, 0); FillVector3D(bottomDV, 0, 0, -1); } } else { if(rotated==false) { FillVector3D(origin, zPosition, width, 0); FillVector3D(rightDV, 0, -1, 0); FillVector3D(bottomDV, 0, 0, 1); } else { FillVector3D(origin, zPosition, 0, height); FillVector3D(rightDV, 0, 1, 0); FillVector3D(bottomDV, 0, 0, -1); } } normalDirection = 0; break; default: itkExceptionMacro("unknown PlaneOrientation"); } if ( transform != NULL ) { origin = transform->TransformPoint( origin ); rightDV = transform->TransformVector( rightDV ); bottomDV = transform->TransformVector( bottomDV ); } ScalarType bounds[6]= { 0, width, 0, height, 0, 1 }; this->SetBounds( bounds ); if ( transform == NULL ) { this->SetMatrixByVectors( rightDV, bottomDV ); } else { this->SetMatrixByVectors( rightDV, bottomDV, transform->GetMatrix().GetVnlMatrix() .get_column(normalDirection).magnitude() ); } this->SetOrigin(origin); } void PlaneGeometry::InitializeStandardPlane( const Geometry3D *geometry3D, PlaneOrientation planeorientation, ScalarType zPosition, bool frontside, bool rotated ) { this->SetReferenceGeometry( const_cast< Geometry3D * >( geometry3D ) ); ScalarType width, height; const BoundingBox::BoundsArrayType& boundsarray = geometry3D->GetBoundingBox()->GetBounds(); Vector3D originVector; FillVector3D(originVector, boundsarray[0], boundsarray[2], boundsarray[4]); + MITK_INFO<GetImageGeometry()) { FillVector3D( originVector, originVector[0] - 0.5, originVector[1] - 0.5, originVector[2] - 0.5 ); } switch(planeorientation) { case Transversal: width = geometry3D->GetExtent(0); height = geometry3D->GetExtent(1); break; case Frontal: width = geometry3D->GetExtent(0); height = geometry3D->GetExtent(2); break; case Sagittal: width = geometry3D->GetExtent(1); height = geometry3D->GetExtent(2); break; default: itkExceptionMacro("unknown PlaneOrientation"); } InitializeStandardPlane( width, height, geometry3D->GetIndexToWorldTransform(), planeorientation, zPosition, frontside, rotated ); ScalarType bounds[6]= { 0, width, 0, height, 0, 1 }; this->SetBounds( bounds ); Point3D origin; originVector = geometry3D->GetIndexToWorldTransform() ->TransformVector( originVector ); origin = GetOrigin() + originVector; SetOrigin(origin); } void PlaneGeometry::InitializeStandardPlane( const Geometry3D *geometry3D, bool top, PlaneOrientation planeorientation, bool frontside, bool rotated ) { ScalarType zPosition; switch(planeorientation) { case Transversal: zPosition = (top ? 0.5 : geometry3D->GetExtent(2)-1+0.5); break; case Frontal: zPosition = (top ? 0.5 : geometry3D->GetExtent(1)-1+0.5); break; case Sagittal: zPosition = (top ? 0.5 : geometry3D->GetExtent(0)-1+0.5); break; default: itkExceptionMacro("unknown PlaneOrientation"); } InitializeStandardPlane( geometry3D, planeorientation, zPosition, frontside, rotated ); } void PlaneGeometry::InitializeStandardPlane( const Vector3D &rightVector, const Vector3D &downVector, const Vector3D *spacing ) { InitializeStandardPlane( rightVector.Get_vnl_vector(), downVector.Get_vnl_vector(), spacing ); } void PlaneGeometry::InitializeStandardPlane( const VnlVector& rightVector, const VnlVector &downVector, const Vector3D *spacing ) { ScalarType width = rightVector.magnitude(); ScalarType height = downVector.magnitude(); InitializeStandardPlane( width, height, rightVector, downVector, spacing ); } void PlaneGeometry::InitializeStandardPlane( mitk::ScalarType width, ScalarType height, const Vector3D &rightVector, const Vector3D &downVector, const Vector3D *spacing ) { InitializeStandardPlane( width, height, rightVector.Get_vnl_vector(), downVector.Get_vnl_vector(), spacing ); } void PlaneGeometry::InitializeStandardPlane( mitk::ScalarType width, ScalarType height, const VnlVector &rightVector, const VnlVector &downVector, const Vector3D *spacing ) { assert(width > 0); assert(height > 0); VnlVector rightDV = rightVector; rightDV.normalize(); VnlVector downDV = downVector; downDV.normalize(); VnlVector normal = vnl_cross_3d(rightVector, downVector); normal.normalize(); if(spacing!=NULL) { rightDV *= (*spacing)[0]; downDV *= (*spacing)[1]; normal *= (*spacing)[2]; } AffineTransform3D::Pointer transform = AffineTransform3D::New(); Matrix3D matrix; matrix.GetVnlMatrix().set_column(0, rightDV); matrix.GetVnlMatrix().set_column(1, downDV); matrix.GetVnlMatrix().set_column(2, normal); transform->SetMatrix(matrix); transform->SetOffset(m_IndexToWorldTransform->GetOffset()); ScalarType bounds[6] = { 0, width, 0, height, 0, 1 }; this->SetBounds( bounds ); this->SetIndexToWorldTransform( transform ); } void PlaneGeometry::InitializePlane( const Point3D &origin, const Vector3D &normal ) { VnlVector rightVectorVnl(3), downVectorVnl; if( Equal( normal[1], 0.0f ) == false ) { FillVector3D( rightVectorVnl, 1.0f, -normal[0]/normal[1], 0.0f ); rightVectorVnl.normalize(); } else { FillVector3D( rightVectorVnl, 0.0f, 1.0f, 0.0f ); } downVectorVnl = vnl_cross_3d( normal.Get_vnl_vector(), rightVectorVnl ); downVectorVnl.normalize(); InitializeStandardPlane( rightVectorVnl, downVectorVnl ); SetOrigin(origin); } void PlaneGeometry::SetMatrixByVectors( const VnlVector &rightVector, const VnlVector &downVector, ScalarType thickness ) { VnlVector normal = vnl_cross_3d(rightVector, downVector); normal.normalize(); normal *= thickness; AffineTransform3D::Pointer transform = AffineTransform3D::New(); Matrix3D matrix; matrix.GetVnlMatrix().set_column(0, rightVector); matrix.GetVnlMatrix().set_column(1, downVector); matrix.GetVnlMatrix().set_column(2, normal); transform->SetMatrix(matrix); transform->SetOffset(m_IndexToWorldTransform->GetOffset()); SetIndexToWorldTransform(transform); } Vector3D PlaneGeometry::GetNormal() const { Vector3D frontToBack; frontToBack.Set_vnl_vector( m_IndexToWorldTransform ->GetMatrix().GetVnlMatrix().get_column(2) ); return frontToBack; } VnlVector PlaneGeometry::GetNormalVnl() const { return m_IndexToWorldTransform ->GetMatrix().GetVnlMatrix().get_column(2); } ScalarType PlaneGeometry::DistanceFromPlane( const Point3D &pt3d_mm ) const { return fabs(SignedDistance( pt3d_mm )); } ScalarType PlaneGeometry::SignedDistance( const Point3D &pt3d_mm ) const { return SignedDistanceFromPlane(pt3d_mm); } bool PlaneGeometry::IsAbove( const Point3D &pt3d_mm ) const { return SignedDistanceFromPlane(pt3d_mm) > 0; } bool PlaneGeometry::IntersectionLine( const PlaneGeometry* plane, Line3D& crossline ) const { Vector3D normal = this->GetNormal(); normal.Normalize(); Vector3D planeNormal = plane->GetNormal(); planeNormal.Normalize(); Vector3D direction = itk::CrossProduct( normal, planeNormal ); if ( direction.GetSquaredNorm() < eps ) return false; crossline.SetDirection( direction ); double N1dN2 = normal * planeNormal; double determinant = 1.0 - N1dN2 * N1dN2; Vector3D origin = this->GetOrigin().GetVectorFromOrigin(); Vector3D planeOrigin = plane->GetOrigin().GetVectorFromOrigin(); double d1 = normal * origin; double d2 = planeNormal * planeOrigin; double c1 = ( d1 - d2 * N1dN2 ) / determinant; double c2 = ( d2 - d1 * N1dN2 ) / determinant; Vector3D p = normal * c1 + planeNormal * c2; crossline.GetPoint().Get_vnl_vector() = p.Get_vnl_vector(); return true; } unsigned int PlaneGeometry::IntersectWithPlane2D( const PlaneGeometry* plane, Point2D& lineFrom, Point2D &lineTo ) const { Line3D crossline; if ( this->IntersectionLine( plane, crossline ) == false ) return 0; Point2D point2; Vector2D direction2; this->Map( crossline.GetPoint(), point2 ); this->Map( crossline.GetPoint(), crossline.GetDirection(), direction2 ); return Line3D::RectangleLineIntersection( 0, 0, GetExtentInMM(0), GetExtentInMM(1), point2, direction2, lineFrom, lineTo ); } double PlaneGeometry::Angle( const PlaneGeometry *plane ) const { return angle(plane->GetMatrixColumn(2), GetMatrixColumn(2)); } double PlaneGeometry::Angle( const Line3D &line ) const { return vnl_math::pi_over_2 - angle( line.GetDirection().Get_vnl_vector(), GetMatrixColumn(2) ); } bool PlaneGeometry::IntersectionPoint( const Line3D &line, Point3D &intersectionPoint ) const { Vector3D planeNormal = this->GetNormal(); planeNormal.Normalize(); Vector3D lineDirection = line.GetDirection(); lineDirection.Normalize(); double t = planeNormal * lineDirection; if ( fabs( t ) < eps ) { return false; } Vector3D diff; diff = this->GetOrigin() - line.GetPoint(); t = ( planeNormal * diff ) / t; intersectionPoint = line.GetPoint() + lineDirection * t; return true; } bool PlaneGeometry::IntersectionPointParam( const Line3D &line, double &t ) const { Vector3D planeNormal = this->GetNormal(); Vector3D lineDirection = line.GetDirection(); t = planeNormal * lineDirection; if ( fabs( t ) < eps ) { return false; } Vector3D diff; diff = this->GetOrigin() - line.GetPoint(); t = ( planeNormal * diff ) / t; return true; } bool PlaneGeometry::IsParallel( const PlaneGeometry *plane ) const { return ( (Angle(plane) < 10.0 * mitk::sqrteps ) || ( Angle(plane) > ( vnl_math::pi - 10.0 * sqrteps ) ) ) ; } bool PlaneGeometry::IsOnPlane( const Point3D &point ) const { return Distance(point) < eps; } bool PlaneGeometry::IsOnPlane( const Line3D &line ) const { return ( (Distance( line.GetPoint() ) < eps) && (Distance( line.GetPoint2() ) < eps) ); } bool PlaneGeometry::IsOnPlane( const PlaneGeometry *plane ) const { return ( IsParallel( plane ) && (Distance( plane->GetOrigin() ) < eps) ); } Point3D PlaneGeometry::ProjectPointOntoPlane( const Point3D& pt ) const { ScalarType len = this->GetNormalVnl().two_norm(); return pt - this->GetNormal() * this->SignedDistanceFromPlane( pt ) / len; } AffineGeometryFrame3D::Pointer PlaneGeometry::Clone() const { Self::Pointer newGeometry = new PlaneGeometry(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } void PlaneGeometry::ExecuteOperation( Operation *operation ) { vtkTransform *transform = vtkTransform::New(); transform->SetMatrix( m_VtkMatrix ); switch ( operation->GetOperationType() ) { case OpORIENT: { mitk::PlaneOperation *planeOp = dynamic_cast< mitk::PlaneOperation * >( operation ); if ( planeOp == NULL ) { return; } Point3D center = planeOp->GetPoint(); Vector3D orientationVector = planeOp->GetNormal(); Vector3D defaultVector; FillVector3D( defaultVector, 0.0, 0.0, 1.0 ); Vector3D rotationAxis = itk::CrossProduct( orientationVector, defaultVector ); //vtkFloatingPointType rotationAngle = acos( orientationVector[2] / orientationVector.GetNorm() ); vtkFloatingPointType rotationAngle = atan2( (double) rotationAxis.GetNorm(), (double) (orientationVector * defaultVector) ); rotationAngle *= 180.0 / vnl_math::pi; transform->PostMultiply(); transform->Identity(); transform->Translate( center[0], center[1], center[2] ); transform->RotateWXYZ( rotationAngle, rotationAxis[0], rotationAxis[1], rotationAxis[2] ); transform->Translate( -center[0], -center[1], -center[2] ); break; } default: Superclass::ExecuteOperation( operation ); transform->Delete(); return; } m_VtkMatrix->DeepCopy(transform->GetMatrix()); this->TransferVtkToItkTransform(); this->Modified(); transform->Delete(); } void PlaneGeometry::PrintSelf( std::ostream& os, itk::Indent indent ) const { Superclass::PrintSelf(os,indent); os << indent << " Normal: " << GetNormal() << std::endl; } } // namespace diff --git a/Modules/MitkExt/Interactions/mitkRegionGrowingTool.cpp b/Modules/MitkExt/Interactions/mitkRegionGrowingTool.cpp index 02afda10ef..1ad71f2e43 100644 --- a/Modules/MitkExt/Interactions/mitkRegionGrowingTool.cpp +++ b/Modules/MitkExt/Interactions/mitkRegionGrowingTool.cpp @@ -1,712 +1,712 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision$ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkRegionGrowingTool.h" #include "mitkToolManager.h" #include "mitkOverwriteSliceImageFilter.h" #include "mitkImageDataItem.h" #include "mitkBaseRenderer.h" #include "mitkRenderingManager.h" #include "mitkApplicationCursor.h" #include "ipSegmentation.h" #include "mitkRegionGrowingTool.xpm" #include "mitkOverwriteDirectedPlaneImageFilter.h" #include "mitkExtractDirectedPlaneImageFilterNew.h" namespace mitk { MITK_TOOL_MACRO(MitkExt_EXPORT, RegionGrowingTool, "Region growing tool"); } #define ROUND(a) ((a)>0 ? (int)((a)+0.5) : -(int)(0.5-(a))) mitk::RegionGrowingTool::RegionGrowingTool() :FeedbackContourTool("PressMoveRelease"), m_LowerThreshold(200), m_UpperThreshold(200), m_InitialLowerThreshold(200), m_InitialUpperThreshold(200), m_ScreenYDifference(0), m_OriginalPicSlice(NULL), m_SeedPointMemoryOffset(0), m_VisibleWindow(0), m_DefaultWindow(0), m_MouseDistanceScaleFactor(3.0), m_LastWorkingSeed(-1), m_FillFeedbackContour(true) { } mitk::RegionGrowingTool::~RegionGrowingTool() { } const char** mitk::RegionGrowingTool::GetXPM() const { return mitkRegionGrowingTool_xpm; } const char* mitk::RegionGrowingTool::GetName() const { return "Region Growing"; } void mitk::RegionGrowingTool::Activated() { Superclass::Activated(); } void mitk::RegionGrowingTool::Deactivated() { Superclass::Deactivated(); } /** 1 Determine which slice is clicked into 2 Determine if the user clicked inside or outside of the segmentation 3 Depending on the pixel value under the mouse click position, two different things happen: (separated out into OnMousePressedInside and OnMousePressedOutside) 3.1 Create a skeletonization of the segmentation and try to find a nice cut 3.1.1 Call a ipSegmentation algorithm to create a nice cut 3.1.2 Set the result of this algorithm as the feedback contour 3.2 Initialize region growing 3.2.1 Determine memory offset inside the original image 3.2.2 Determine initial region growing parameters from the level window settings of the image 3.2.3 Perform a region growing (which generates a new feedback contour) */ bool mitk::RegionGrowingTool::OnMousePressed (Action* action, const StateEvent* stateEvent) { //ToolLogger::SetVerboseness(3); MITK_INFO << "OnMousePressed" << std::endl; if (FeedbackContourTool::OnMousePressed( action, stateEvent )) { MITK_INFO << "OnMousePressed: FeedbackContourTool says ok" << std::endl; // 1. Find out which slice the user clicked, find out which slice of the toolmanager's reference and working image corresponds to that const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (positionEvent) { MITK_INFO << "OnMousePressed: got positionEvent" << std::endl; m_ReferenceSlice = FeedbackContourTool::GetAffectedReferenceSlice( positionEvent ); m_WorkingSlice = FeedbackContourTool::GetAffectedWorkingSlice( positionEvent ); if ( m_WorkingSlice.IsNotNull() ) // can't do anything without the segmentation { MITK_INFO << "OnMousePressed: got working slice" << std::endl; // 2. Determine if the user clicked inside or outside of the segmentation const Geometry3D* workingSliceGeometry = m_WorkingSlice->GetGeometry(); Point3D mprojectedPointIn2D; workingSliceGeometry->WorldToIndex( positionEvent->GetWorldPosition(), mprojectedPointIn2D); itk::Index<2> projectedPointInWorkingSlice2D; projectedPointInWorkingSlice2D[0] = static_cast( mprojectedPointIn2D[0] - 0.5 ); projectedPointInWorkingSlice2D[1] = static_cast( mprojectedPointIn2D[1] - 0.5 ); if ( workingSliceGeometry->IsIndexInside( projectedPointInWorkingSlice2D ) ) { MITK_INFO << "OnMousePressed: point " << positionEvent->GetWorldPosition() << " (index coordinates " << projectedPointInWorkingSlice2D << ") IS in working slice" << std::endl; // Convert to ipMITKSegmentationTYPE (because getting pixels relys on that data type) itk::Image< ipMITKSegmentationTYPE, 2 >::Pointer correctPixelTypeImage; CastToItkImage( m_WorkingSlice, correctPixelTypeImage ); assert (correctPixelTypeImage.IsNotNull() ); // possible bug in CastToItkImage ? // direction maxtrix is wrong/broken/not working after CastToItkImage, leading to a failed assertion in // mitk/Core/DataStructures/mitkSlicedGeometry3D.cpp, 479: // virtual void mitk::SlicedGeometry3D::SetSpacing(const mitk::Vector3D&): Assertion `aSpacing[0]>0 && aSpacing[1]>0 && aSpacing[2]>0' failed // solution here: we overwrite it with an unity matrix itk::Image< ipMITKSegmentationTYPE, 2 >::DirectionType imageDirection; imageDirection.SetIdentity(); correctPixelTypeImage->SetDirection(imageDirection); Image::Pointer temporarySlice = Image::New(); // temporarySlice = ImportItkImage( correctPixelTypeImage ); CastToMitkImage( correctPixelTypeImage, temporarySlice ); mitkIpPicDescriptor* workingPicSlice = temporarySlice->GetSliceData()->GetPicDescriptor(); int initialWorkingOffset = projectedPointInWorkingSlice2D[1] * workingPicSlice->n[0] + projectedPointInWorkingSlice2D[0]; if ( initialWorkingOffset < static_cast( workingPicSlice->n[0] * workingPicSlice->n[1] ) && initialWorkingOffset >= 0 ) { // 3. determine the pixel value under the last click bool inside = static_cast(workingPicSlice->data)[initialWorkingOffset] != 0; m_PaintingPixelValue = inside ? 0 : 1; // if inside, we want to remove a part, otherwise we want to add something if ( m_LastWorkingSeed >= static_cast( workingPicSlice->n[0] * workingPicSlice->n[1] ) || m_LastWorkingSeed < 0 ) { inside = false; } if ( m_ReferenceSlice.IsNotNull() ) { MITK_INFO << "OnMousePressed: got reference slice" << std::endl; m_OriginalPicSlice = m_ReferenceSlice->GetSliceData()->GetPicDescriptor(); // 3.1. Switch depending on the pixel value if (inside) { OnMousePressedInside(action, stateEvent, workingPicSlice, initialWorkingOffset); } else { OnMousePressedOutside(action, stateEvent); } } } } } } } MITK_INFO << "end OnMousePressed" << std::endl; return true; } /** 3.1 Create a skeletonization of the segmentation and try to find a nice cut 3.1.1 Call a ipSegmentation algorithm to create a nice cut 3.1.2 Set the result of this algorithm as the feedback contour */ bool mitk::RegionGrowingTool::OnMousePressedInside(Action* itkNotUsed( action ), const StateEvent* stateEvent, mitkIpPicDescriptor* workingPicSlice, int initialWorkingOffset) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); // checked in OnMousePressed // 3.1.1. Create a skeletonization of the segmentation and try to find a nice cut // apply the skeletonization-and-cut algorithm // generate contour to remove // set m_ReferenceSlice = NULL so nothing will happen during mouse move // remember to fill the contour with 0 in mouserelease mitkIpPicDescriptor* segmentationHistory = ipMITKSegmentationCreateGrowerHistory( workingPicSlice, m_LastWorkingSeed, NULL ); // free again if (segmentationHistory) { tCutResult cutContour = ipMITKSegmentationGetCutPoints( workingPicSlice, segmentationHistory, initialWorkingOffset ); // tCutResult is a ipSegmentation type mitkIpPicFree( segmentationHistory ); if (cutContour.cutIt) { // 3.1.2 copy point from float* to mitk::Contour Contour::Pointer contourInImageIndexCoordinates = Contour::New(); contourInImageIndexCoordinates->Initialize(); Point3D newPoint; for (int index = 0; index < cutContour.deleteSize; ++index) { newPoint[0] = cutContour.deleteCurve[ 2 * index + 0 ]; newPoint[1] = cutContour.deleteCurve[ 2 * index + 1 ]; newPoint[2] = 0.0; contourInImageIndexCoordinates->AddVertex( newPoint - 0.5 ); } free(cutContour.traceline); free(cutContour.deleteCurve); // perhaps visualize this for fun? free(cutContour.onGradient); Contour::Pointer contourInWorldCoordinates = FeedbackContourTool::BackProjectContourFrom2DSlice( m_WorkingSlice->GetGeometry(), contourInImageIndexCoordinates, true ); // true: sub 0.5 for ipSegmentation correction FeedbackContourTool::SetFeedbackContour( *contourInWorldCoordinates ); FeedbackContourTool::SetFeedbackContourVisible(true); mitk::RenderingManager::GetInstance()->RequestUpdate( positionEvent->GetSender()->GetRenderWindow() ); m_FillFeedbackContour = true; } else { m_FillFeedbackContour = false; } } else { m_FillFeedbackContour = false; } m_ReferenceSlice = NULL; return true; } /** 3.2 Initialize region growing 3.2.1 Determine memory offset inside the original image 3.2.2 Determine initial region growing parameters from the level window settings of the image 3.2.3 Perform a region growing (which generates a new feedback contour) */ bool mitk::RegionGrowingTool::OnMousePressedOutside(Action* itkNotUsed( action ), const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); // checked in OnMousePressed // 3.2 If we have a reference image, then perform an initial region growing, considering the reference image's level window // if click was outside the image, don't continue const Geometry3D* sliceGeometry = m_ReferenceSlice->GetGeometry(); Point3D mprojectedPointIn2D; sliceGeometry->WorldToIndex( positionEvent->GetWorldPosition(), mprojectedPointIn2D ); itk::Index<2> projectedPointIn2D; projectedPointIn2D[0] = static_cast( mprojectedPointIn2D[0] - 0.5 ); projectedPointIn2D[1] = static_cast( mprojectedPointIn2D[1] - 0.5 ); - if ( sliceGeometry->IsIndexInside( mprojectedPointIn2D ) ) + if ( sliceGeometry->IsIndexInside( projectedPointIn2D ) ) { MITK_INFO << "OnMousePressed: point " << positionEvent->GetWorldPosition() << " (index coordinates " << mprojectedPointIn2D << ") IS in reference slice" << std::endl; // 3.2.1 Remember Y cursor position and initial seed point //m_ScreenYPositionAtStart = static_cast(positionEvent->GetDisplayPosition()[1]); m_LastScreenPosition = ApplicationCursor::GetInstance()->GetCursorPosition(); m_ScreenYDifference = 0; m_SeedPointMemoryOffset = projectedPointIn2D[1] * m_OriginalPicSlice->n[0] + projectedPointIn2D[0]; m_LastWorkingSeed = m_SeedPointMemoryOffset; // remember for skeletonization if ( m_SeedPointMemoryOffset < static_cast( m_OriginalPicSlice->n[0] * m_OriginalPicSlice->n[1] ) && m_SeedPointMemoryOffset >= 0 ) { // 3.2.2 Get level window from reference DataNode // Use some logic to determine initial gray value bounds LevelWindow lw(0, 500); m_ToolManager->GetReferenceData(0)->GetLevelWindow(lw); // will fill lw if levelwindow property is present, otherwise won't touch it. m_VisibleWindow = lw.GetWindow(); // necessary for limiting the upper and lower threshold to the maximum gray values m_DefaultWindow = lw.GetDefaultWindow(); static bool initializedAlready = false; // just evaluated once if (!initializedAlready) { m_InitialLowerThreshold = static_cast(m_VisibleWindow / 10.0); // 20% of the visible gray values m_InitialUpperThreshold = static_cast(m_VisibleWindow / 10.0); initializedAlready = true; } m_LowerThreshold = m_InitialLowerThreshold; m_UpperThreshold = m_InitialUpperThreshold; DisplayGeometry* displayGeometry = positionEvent->GetSender()->GetDisplayGeometry(); if (displayGeometry) { m_MouseDistanceScaleFactor = m_VisibleWindow / ( 3.0 * displayGeometry->GetDisplayHeight() ); } // 3.2.3. Actually perform region growing mitkIpPicDescriptor* result = PerformRegionGrowingAndUpdateContour(); ipMITKSegmentationFree( result); // display the contour FeedbackContourTool::SetFeedbackContourVisible(true); mitk::RenderingManager::GetInstance()->RequestUpdate(positionEvent->GetSender()->GetRenderWindow()); m_FillFeedbackContour = true; } } return true; } /** If in region growing mode (m_ReferenceSlice != NULL), then 1. Calculate the new thresholds from mouse position (relative to first position) 2. Perform a new region growing and update the feedback contour */ bool mitk::RegionGrowingTool::OnMouseMoved (Action* action, const StateEvent* stateEvent) { if (FeedbackContourTool::OnMouseMoved( action, stateEvent )) { if ( m_ReferenceSlice.IsNotNull() && m_OriginalPicSlice ) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (positionEvent) { ApplicationCursor* cursor = ApplicationCursor::GetInstance(); if (!cursor) return false; m_ScreenYDifference += cursor->GetCursorPosition()[1] - m_LastScreenPosition[1]; cursor->SetCursorPosition( m_LastScreenPosition ); m_LowerThreshold = m_InitialLowerThreshold + static_cast( m_ScreenYDifference * m_MouseDistanceScaleFactor ); if (m_LowerThreshold < 1) m_LowerThreshold = 1; if (m_LowerThreshold > m_VisibleWindow / 2) m_LowerThreshold = m_VisibleWindow / 2; m_UpperThreshold = m_InitialUpperThreshold + static_cast( m_ScreenYDifference * m_MouseDistanceScaleFactor ); if (m_UpperThreshold < 1) m_UpperThreshold = 1; if (m_UpperThreshold > m_VisibleWindow / 2) m_UpperThreshold = m_VisibleWindow / 2; //MITK_INFO << "new interval: l " << m_LowerThreshold << " u " << m_UpperThreshold << std::endl; // 2. Perform region growing again and show the result mitkIpPicDescriptor* result = PerformRegionGrowingAndUpdateContour(); ipMITKSegmentationFree( result ); // 3. Update the contour mitk::RenderingManager::GetInstance()->ForceImmediateUpdate(positionEvent->GetSender()->GetRenderWindow()); } } } return true; } /** If the feedback contour should be filled, then it is done here. (Contour is NOT filled, when skeletonization is done but no nice cut was found) */ bool mitk::RegionGrowingTool::OnMouseReleased(Action* action, const StateEvent* stateEvent) { if (FeedbackContourTool::OnMouseReleased( action, stateEvent )) { // 1. If we have a working slice, use the contour to fill a new piece on segmentation on it (or erase a piece that was selected by ipMITKSegmentationGetCutPoints) if ( m_WorkingSlice.IsNotNull() && m_OriginalPicSlice ) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (positionEvent) { // remember parameters for next time m_InitialLowerThreshold = m_LowerThreshold; m_InitialUpperThreshold = m_UpperThreshold; if (m_FillFeedbackContour) { // 3. use contour to fill a region in our working slice Contour* feedbackContour( FeedbackContourTool::GetFeedbackContour() ); if (feedbackContour) { Contour::Pointer projectedContour = FeedbackContourTool::ProjectContourTo2DSlice( m_WorkingSlice, feedbackContour, false, false ); // false: don't add any 0.5 // false: don't constrain the contour to the image's inside if (projectedContour.IsNotNull()) { FeedbackContourTool::FillContourInSlice( projectedContour, m_WorkingSlice, m_PaintingPixelValue ); // 4. write working slice back into image volume int affectedDimension( -1 ); int affectedSlice( -1 ); const PlaneGeometry* planeGeometry( dynamic_cast (positionEvent->GetSender()->GetCurrentWorldGeometry2D() ) ); FeedbackContourTool::DetermineAffectedImageSlice( dynamic_cast( m_ToolManager->GetWorkingData(0)->GetData() ), planeGeometry, affectedDimension, affectedSlice ); MITK_INFO << "OnMouseReleased: writing back to dimension " << affectedDimension << ", slice " << affectedSlice << " in working image" << std::endl; //If dazu gemacht Image::Pointer workingImage = dynamic_cast( m_ToolManager->GetWorkingData(0)->GetData() ); if ( affectedDimension != -1 ) { OverwriteSliceImageFilter::Pointer slicewriter = OverwriteSliceImageFilter::New(); slicewriter->SetInput( workingImage ); slicewriter->SetCreateUndoInformation( true ); slicewriter->SetSliceImage( m_WorkingSlice ); slicewriter->SetSliceDimension( affectedDimension ); slicewriter->SetSliceIndex( affectedSlice ); slicewriter->SetTimeStep( positionEvent->GetSender()->GetTimeStep( workingImage ) ); slicewriter->Update(); if ( m_RememberContourPositions ) { this->AddContourmarker(positionEvent); } } else { OverwriteDirectedPlaneImageFilter::Pointer slicewriter = OverwriteDirectedPlaneImageFilter::New(); slicewriter->SetInput( workingImage ); slicewriter->SetCreateUndoInformation( false ); slicewriter->SetSliceImage( m_WorkingSlice ); slicewriter->SetPlaneGeometry3D( m_WorkingSlice->GetGeometry() ); slicewriter->SetTimeStep( positionEvent->GetSender()->GetTimeStep( workingImage ) ); slicewriter->Update(); if ( m_RememberContourPositions ) { this->AddContourmarker(positionEvent); } } } } } FeedbackContourTool::SetFeedbackContourVisible(false); mitk::RenderingManager::GetInstance()->RequestUpdate( positionEvent->GetSender()->GetRenderWindow() ); } } } m_ReferenceSlice = NULL; // don't leak m_WorkingSlice = NULL; m_OriginalPicSlice = NULL; return true; } /** Uses ipSegmentation algorithms to do the actual region growing. The result (binary image) is first smoothed by a 5x5 circle mask, then its contour is extracted and converted to MITK coordinates. */ mitkIpPicDescriptor* mitk::RegionGrowingTool::PerformRegionGrowingAndUpdateContour() { // 1. m_OriginalPicSlice and m_SeedPointMemoryOffset are set to sensitive values, as well as m_LowerThreshold and m_UpperThreshold assert (m_OriginalPicSlice); if (m_OriginalPicSlice->n[0] != 256 || m_OriginalPicSlice->n[1] != 256) // ??? assert( (m_SeedPointMemoryOffset < static_cast( m_OriginalPicSlice->n[0] * m_OriginalPicSlice->n[1] )) && (m_SeedPointMemoryOffset >= 0) ); // inside the image // 2. ipSegmentation is used to perform region growing float ignored; int oneContourOffset( 0 ); mitkIpPicDescriptor* regionGrowerResult = ipMITKSegmentationGrowRegion4N( m_OriginalPicSlice, m_SeedPointMemoryOffset, // seed point true, // grayvalue interval relative to seed point gray value? m_LowerThreshold, m_UpperThreshold, 0, // continue until done (maxIterations == 0) NULL, // allocate new memory (only this time, on mouse move we'll reuse the old buffer) oneContourOffset, // a pixel that is near the resulting contour ignored // ignored by us ); if (!regionGrowerResult || oneContourOffset == -1) { Contour::Pointer dummyContour = Contour::New(); dummyContour->Initialize(); FeedbackContourTool::SetFeedbackContour( *dummyContour ); if (regionGrowerResult) ipMITKSegmentationFree(regionGrowerResult); return NULL; } // 3. We smooth the result a little to reduce contour complexity bool smoothResult( true ); // currently fixed, perhaps remove else block mitkIpPicDescriptor* smoothedRegionGrowerResult; if (smoothResult) { // Smooth the result (otherwise very detailed contour) smoothedRegionGrowerResult = SmoothIPPicBinaryImage( regionGrowerResult, oneContourOffset ); ipMITKSegmentationFree( regionGrowerResult ); } else { smoothedRegionGrowerResult = regionGrowerResult; } // 4. convert the result of region growing into a mitk::Contour // At this point oneContourOffset could be useless, if smoothing destroyed a thin bridge. In these // cases, we have two or more unconnected segmentation regions, and we don't know, which one is touched by oneContourOffset. // In the bad case, the contour is not the one around our seedpoint, so the result looks very strange to the user. // -> we remove the point where the contour started so far. Then we look from the bottom of the image for the first segmentation pixel // and start another contour extraction from there. This is done, until the seedpoint is inside the contour int numberOfContourPoints( 0 ); int newBufferSize( 0 ); float* contourPoints = ipMITKSegmentationGetContour8N( smoothedRegionGrowerResult, oneContourOffset, numberOfContourPoints, newBufferSize ); // memory allocated with malloc if (contourPoints) { while ( !ipMITKSegmentationIsInsideContour( contourPoints, // contour numberOfContourPoints, // points in contour m_SeedPointMemoryOffset % smoothedRegionGrowerResult->n[0], // test point x m_SeedPointMemoryOffset / smoothedRegionGrowerResult->n[0] // test point y ) ) { // we decide that this cannot be part of the segmentation because the seedpoint is not contained in the contour (fill the 4-neighborhood with 0) ipMITKSegmentationReplaceRegion4N( smoothedRegionGrowerResult, oneContourOffset, 0 ); // move the contour offset to the last row (x position of the seed point) int rowLength = smoothedRegionGrowerResult->n[0]; // number of pixels in a row oneContourOffset = m_SeedPointMemoryOffset % smoothedRegionGrowerResult->n[0] // x of seed point + rowLength*(smoothedRegionGrowerResult->n[1]-1); // y of last row while ( oneContourOffset >=0 && (*(static_cast(smoothedRegionGrowerResult->data) + oneContourOffset) == 0) ) { oneContourOffset -= rowLength; // if pixel at data+oneContourOffset is 0, then move up one row } if ( oneContourOffset < 0 ) { break; // just use the last contour we found } free(contourPoints); // release contour memory contourPoints = ipMITKSegmentationGetContour8N( smoothedRegionGrowerResult, oneContourOffset, numberOfContourPoints, newBufferSize ); // memory allocated with malloc } // copy point from float* to mitk::Contour Contour::Pointer contourInImageIndexCoordinates = Contour::New(); contourInImageIndexCoordinates->Initialize(); Point3D newPoint; for (int index = 0; index < numberOfContourPoints; ++index) { newPoint[0] = contourPoints[ 2 * index + 0 ]; newPoint[1] = contourPoints[ 2 * index + 1 ]; newPoint[2] = 0; contourInImageIndexCoordinates->AddVertex( newPoint - 0.5); } free(contourPoints); Contour::Pointer contourInWorldCoordinates = FeedbackContourTool::BackProjectContourFrom2DSlice( m_ReferenceSlice->GetGeometry(), contourInImageIndexCoordinates, true ); // true: sub 0.5 for ipSegmentation correctio FeedbackContourTool::SetFeedbackContour( *contourInWorldCoordinates ); } // 5. Result HAS TO BE freed by caller, contains the binary region growing result return smoothedRegionGrowerResult; } /** Helper method for SmoothIPPicBinaryImage. Smoothes a given part of and image. \param sourceImage The original binary image. \param dest The smoothed image (will be written without bounds checking). \param contourOfs One offset of the contour. Is updated if a pixel is changed (which might change the contour). \param maskOffsets Memory offsets that describe the smoothing mask. \param maskSize Entries of the mask. \param startOffset First pixel that should be smoothed using this mask. \param endOffset Last pixel that should be smoothed using this mask. */ void mitk::RegionGrowingTool::SmoothIPPicBinaryImageHelperForRows( mitkIpPicDescriptor* sourceImage, mitkIpPicDescriptor* dest, int &contourOfs, int* maskOffsets, int maskSize, int startOffset, int endOffset ) { // work on the very first row ipMITKSegmentationTYPE* current; ipMITKSegmentationTYPE* source = ((ipMITKSegmentationTYPE*)sourceImage->data) + startOffset; // + 1! don't read at start-1 ipMITKSegmentationTYPE* end = ((ipMITKSegmentationTYPE*)dest->data) + endOffset; int ofs = startOffset; int minority = (maskSize - 1) / 2; for (current = ((ipMITKSegmentationTYPE*)dest->data) + startOffset; current minority) { *current = 1; contourOfs = ofs; } else { *current = 0; } ++source; ++ofs; } } /** Smoothes a binary ipPic image with a 5x5 mask. The image borders (some first and last rows) are treated differently. */ mitkIpPicDescriptor* mitk::RegionGrowingTool::SmoothIPPicBinaryImage( mitkIpPicDescriptor* image, int &contourOfs, mitkIpPicDescriptor* dest ) { if (!image) return NULL; // Original code from /trunk/mbi-qm/Qmitk/Qmitk2DSegTools/RegionGrowerTool.cpp (first version by T. Boettger?). Reformatted and documented and restructured. #define MSK_SIZE5x5 21 #define MSK_SIZE3x3 5 #define MSK_SIZE3x1 3 // mask is an array of coordinates that form a rastered circle like this // // OOO // OOOOO // OOOOO // OOOOO // OOO // // int mask5x5[MSK_SIZE5x5][2] = { /******/ {-1,-2}, {0,-2}, {1,-2}, /*****/ {-2,-1}, {-1,-1}, {0,-1}, {1,-1}, {2,-1}, {-2, 0}, {-1, 0}, {0, 0}, {1, 0}, {2, 0}, {-2, 1}, {-1, 1}, {0, 1}, {1, 1}, {2, 1}, /******/ {-1, 2}, {0, 2}, {1, 2} /*****/ }; int mask3x3[MSK_SIZE3x3][2] = { /******/ {0,-1}, /*****/ {-1, 0}, {0, 0}, {1, 0}, /******/ {0, 1} /*****/ }; int mask3x1[MSK_SIZE3x1][2] = { {-1, 0}, {0, 0}, {1, 0} }; // The following lines iterate over all the pixels of a (sliced) image (except the first and last three rows). // For each pixel, all the coordinates around it (according to mask) are evaluated (this means 21 pixels). // If more than 10 of the evaluated pixels are non-zero, then the central pixel is set to 1, else to 0. // This is determining a majority. If there is no clear majority, then the central pixel itself "decides". int maskOffset5x5[MSK_SIZE5x5]; int line = image->n[0]; for (int i=0; in[0]; int spareOut1Rows = 1*image->n[0]; if ( image->n[1] > 0 ) SmoothIPPicBinaryImageHelperForRows( image, dest, contourOfs, maskOffset3x1, MSK_SIZE3x1, 1, dest->n[0] ); if ( image->n[1] > 3 ) SmoothIPPicBinaryImageHelperForRows( image, dest, contourOfs, maskOffset3x3, MSK_SIZE3x3, spareOut1Rows, dest->n[0]*3 ); if ( image->n[1] > 6 ) SmoothIPPicBinaryImageHelperForRows( image, dest, contourOfs, maskOffset5x5, MSK_SIZE5x5, spareOut3Rows, dest->n[0]*dest->n[1] - spareOut3Rows ); if ( image->n[1] > 8 ) SmoothIPPicBinaryImageHelperForRows( image, dest, contourOfs, maskOffset3x3, MSK_SIZE3x3, dest->n[0]*dest->n[1] -spareOut3Rows, dest->n[0]*dest->n[1] - spareOut1Rows ); if ( image->n[1] > 10) SmoothIPPicBinaryImageHelperForRows( image, dest, contourOfs, maskOffset3x1, MSK_SIZE3x1, dest->n[0]*dest->n[1] -spareOut1Rows, dest->n[0]*dest->n[1] - 1 ); // correction for first pixel (sorry for the ugliness) if ( *((ipMITKSegmentationTYPE*)(dest->data)+1) == 1 ) { *((ipMITKSegmentationTYPE*)(dest->data)+0) = 1; } if (dest->n[0] * dest->n[1] > 2) { // correction for last pixel if ( *((ipMITKSegmentationTYPE*)(dest->data)+dest->n[0]*dest->n[1]-2) == 1 ) { *((ipMITKSegmentationTYPE*)(dest->data)+dest->n[0]*dest->n[1]-1) = 1; } } return dest; }