Uploaded image for project: 'Data Management'
  1. Data Management
  2. DM-9675

Soften symmetry operator requirements

    Details

    • Type: Story
    • Status: Done
    • Resolution: Done
    • Fix Version/s: None
    • Component/s: meas_deblender
    • Labels:
      None

      Description

      The current symmetry operator enforces strict symmetry, meaning that all pixels that do not have a symmetric partner in the footprint are penalized. To ease this requirement, which may be necessary if some of the flux for one (or more) of the objects lies outside of the footprint, we allow the user to use a less stringent penalty value.

        Attachments

          Activity

          Hide
          fred3m Fred Moolekamp added a comment -

          For simplicity this work is not being done on a separate branch, but instead will be included in DM-9172. The code for the new symmetry operator appears below:

          def getPeakSymmetry(shape, px, py, fillValue=0):
              """Build the operator to symmetrize a the intensities for a single row
              """
              center = (np.array(shape)-1)/2.0
              # If the peak is centered at the middle of the footprint,
              # make the entire footprint symmetric
              if px==center[1] and py==center[0]:
                  return scipy.sparse.coo_matrix(np.fliplr(np.eye(shape[0]*shape[1])))
           
              # Otherwise, find the bounding box that contains the minimum number of pixels needed to symmetrize
              if py<(shape[0]-1)/2.:
                  ymin = 0
                  ymax = 2*py+1
              elif py>(shape[0]-1)/2.:
                  ymin = 2*py-shape[0]+1
                  ymax = shape[0]
              else:
                  ymin = 0
                  ymax = shape[0]
              if px<(shape[1]-1)/2.:
                  xmin = 0
                  xmax = 2*px+1
              elif px>(shape[1]-1)/2.:
                  xmin = 2*px-shape[1]+1
                  xmax = shape[1]
              else:
                  xmin = 0
                  xmax = shape[1]
           
              fpHeight, fpWidth = shape
              fpSize = fpWidth*fpHeight
              tWidth = xmax-xmin
              tHeight = ymax-ymin
              extraWidth = fpWidth-tWidth
              pixels = (tHeight-1)*fpWidth+tWidth
           
              # This is the block of the matrix that symmetrizes intensities at the peak position
              subOp = np.eye(pixels, pixels)
              for i in range(0,tHeight-1):
                  for j in range(extraWidth):
                      idx = (i+1)*tWidth+(i*extraWidth)+j
                      subOp[idx, idx] = fillValue
              subOp = np.fliplr(subOp)
           
              smin = ymin*fpWidth+xmin
              smax = (ymax-1)*fpWidth+xmax
              if fillValue!=0:
                  symmetryOp = np.identity(fpSize)*fillValue
              else:
                  symmetryOp = np.zeros((fpSize, fpSize))
              symmetryOp[smin:smax,smin:smax] = subOp
           
              # Return a sparse matrix, which greatly speeds up the processing
              return scipy.sparse.coo_matrix(symmetryOp)
           
          def getPeakSymmetryOp(shape, px, py, fillValue=0):
              """Operator to calculate the difference from the symmetric intensities
              """
              symOp = getPeakSymmetry(shape, px, py, fillValue)
              diffOp = scipy.sparse.identity(symOp.shape[0])-symOp
              # In cases where the symmetry operator is very small (eg. a small isolated source)
              # scipy doesn't return a sparse matrix, so we test whether or not the matrix is sparse
              # and if it is, use a sparse matrix that works best with the proximal operators.
              if hasattr(diffOp, "tocoo"):
                  diffOp = diffOp.tocoo()
              return diffOp
          

          Show
          fred3m Fred Moolekamp added a comment - For simplicity this work is not being done on a separate branch, but instead will be included in DM-9172 . The code for the new symmetry operator appears below: def getPeakSymmetry(shape, px, py, fillValue = 0 ): """Build the operator to symmetrize a the intensities for a single row """ center = (np.array(shape) - 1 ) / 2.0 # If the peak is centered at the middle of the footprint, # make the entire footprint symmetric if px = = center[ 1 ] and py = = center[ 0 ]: return scipy.sparse.coo_matrix(np.fliplr(np.eye(shape[ 0 ] * shape[ 1 ])))   # Otherwise, find the bounding box that contains the minimum number of pixels needed to symmetrize if py<(shape[ 0 ] - 1 ) / 2. : ymin = 0 ymax = 2 * py + 1 elif py>(shape[ 0 ] - 1 ) / 2. : ymin = 2 * py - shape[ 0 ] + 1 ymax = shape[ 0 ] else : ymin = 0 ymax = shape[ 0 ] if px<(shape[ 1 ] - 1 ) / 2. : xmin = 0 xmax = 2 * px + 1 elif px>(shape[ 1 ] - 1 ) / 2. : xmin = 2 * px - shape[ 1 ] + 1 xmax = shape[ 1 ] else : xmin = 0 xmax = shape[ 1 ]   fpHeight, fpWidth = shape fpSize = fpWidth * fpHeight tWidth = xmax - xmin tHeight = ymax - ymin extraWidth = fpWidth - tWidth pixels = (tHeight - 1 ) * fpWidth + tWidth   # This is the block of the matrix that symmetrizes intensities at the peak position subOp = np.eye(pixels, pixels) for i in range ( 0 ,tHeight - 1 ): for j in range (extraWidth): idx = (i + 1 ) * tWidth + (i * extraWidth) + j subOp[idx, idx] = fillValue subOp = np.fliplr(subOp)   smin = ymin * fpWidth + xmin smax = (ymax - 1 ) * fpWidth + xmax if fillValue! = 0 : symmetryOp = np.identity(fpSize) * fillValue else : symmetryOp = np.zeros((fpSize, fpSize)) symmetryOp[smin:smax,smin:smax] = subOp   # Return a sparse matrix, which greatly speeds up the processing return scipy.sparse.coo_matrix(symmetryOp)   def getPeakSymmetryOp(shape, px, py, fillValue = 0 ): """Operator to calculate the difference from the symmetric intensities """ symOp = getPeakSymmetry(shape, px, py, fillValue) diffOp = scipy.sparse.identity(symOp.shape[ 0 ]) - symOp # In cases where the symmetry operator is very small (eg. a small isolated source) # scipy doesn't return a sparse matrix, so we test whether or not the matrix is sparse # and if it is, use a sparse matrix that works best with the proximal operators. if hasattr (diffOp, "tocoo" ): diffOp = diffOp.tocoo() return diffOp
          Hide
          fred3m Fred Moolekamp added a comment -

          Hi Peter,
          Would you mind taking a look at the new symmetry operator and testing it out?

          It requires you to add a new optional parameter fillValue=0 to nmf_deblender, and call getPeakSymmetryOp using C = getPeakSymmetryOp((N,M), px, py, fillValue). fillValue is the inverse of the penalty to use, so fillValue=0 is the current implementation that strictly enforces symmetry, while fillValue=1 doesn't penalize pixels outside of the symmetric patch at all.

          Let me know if you have any problems or questions.

          -Fred

          Show
          fred3m Fred Moolekamp added a comment - Hi Peter, Would you mind taking a look at the new symmetry operator and testing it out? It requires you to add a new optional parameter fillValue=0 to nmf_deblender , and call getPeakSymmetryOp using C = getPeakSymmetryOp((N,M), px, py, fillValue) . fillValue is the inverse of the penalty to use, so fillValue=0 is the current implementation that strictly enforces symmetry, while fillValue=1 doesn't penalize pixels outside of the symmetric patch at all. Let me know if you have any problems or questions. -Fred
          Hide
          pmelchior Peter Melchior added a comment -

          Works as expected. Interesting to note: fillValue can be some float, but it's only effective in a narrow range around 1.

          Show
          pmelchior Peter Melchior added a comment - Works as expected. Interesting to note: fillValue can be some float, but it's only effective in a narrow range around 1.

            People

            • Assignee:
              fred3m Fred Moolekamp
              Reporter:
              fred3m Fred Moolekamp
              Reviewers:
              Peter Melchior
              Watchers:
              Fred Moolekamp, Peter Melchior
            • Votes:
              0 Vote for this issue
              Watchers:
              2 Start watching this issue

              Dates

              • Created:
                Updated:
                Resolved:

                Summary Panel